Geochemical evolution of groundwater in the Culebra dolomite near the Waste Isolation Pilot Plant, southeastern New Mexico, USA

Geochimica et Cosmochimica Acta, Vol. 58, No. 10, pp. 2299-2323, 1994

Copyright0 1994ElsevierScienceLtd Printedin the USA. All rights reserved

Pergamon

0016-7037/94 $6.00 + .OO

0016-7037(94)EOO39-N

Geochemical evolution of groundwater in the Culebra Dolomite near the Waste Isolation Pilot Plant, southeastern New Mexico, USA* MALCOLMD. SIEGEL’ and SCOTT ANDERHOLM’ ‘Geoscience and Geotechnology Center, Sandia National Laboratories, Albuquerque, NM 87 185, USA 2Water Resources Division, US Geological Survey, Albuquerque, NM 87110, USA (Received April 14, 1993; accepted in revisedform December 29, I993 )

Ah&act-The Culebra Dolomite Member of the Rustler Formation, a thin ( 10 m) fractured dolomite aquifer, lies approximately 450 m above the repository horizon of the Waste Isolation Pilot Plant ( WIPP) in southeastern New Mexico, USA. Salinities of water in the Culebra range roughly from 10,000 to 200,000 mg/L within the WIPP site. A proposed model for the post-Pleistocene hydrochemical evolution of the Culebra tentatively identifies the major sources and sinks for many of the groundwater solutes. Reaction-path simulations with the PHRQPITZ code suggest that the Culebra dolomite is a partial chemical equilibrium system whose composition is controlled by an irreversible process (dissolution of evaporites) and equilibrium with gypsum and calcite. Net geochemical reactions along postulated modem flow paths, calculated with the NETPATH code, include dissolution of halite, carbonate and evaporite salts, and ion exchange. R-mode principal component analysis revealed correlations among the concentrations of Si, Mg, pH, Li, and B that are consistent with several clay-water reactions. The results of the geochemical calculations and mineralogical data are consistent with the following hydrochemical model: ( 1) solutes are added to the Culebra by dissolution of evaporite minerals; (2) the solubilities of gypsum and calcite increase as the salinity increases; these minerals dissolve as chemical equilibrium is maintained between them and the groundwater; (3) equilibrium is not maintained between the waters and dolomite; sufficient Mg is added to the waters by dissolution of accessory camallite or polyhalite such that the degree of dolomite supersaturation increases with ionic strength; and (4) clays within the fractures and rock matrix exert some control on the distribution of Li, B, Mg, and Si via sorption, ion exchange, and dissolution. INTRODUCTION THE WASTE ISOLATIONPILOT PLANT ( WIPP) is being designed for the geologic disposal of defense-generated transuranic (TRU) wastes within the Salado Formation, a thick Ochoan (Permian) evaporite deposit in the northern Delaware Basin of southeastern New Mexico, USA. The Culebra dolomite is a thin (10 m), permeable, laterally continuous unit that lies approximately 450 m above the waste repository horizon and is the most likely pathway for long-term migration of radionuclides to the accessible environment. An understanding of the geochemical evolution of waters in the Culebra Dolomite Member of the overlying Rustler Formation is important to the evaluation of the ability of the bedded evaporite environment to contain waste radionuclides for long periods of time. The hydrogeochemical system of the Culebra and adjacent units in the Rustler Formation has been the subject of intense study since 1976 (see, e.g., MERCER, 1983; GONZALEZ,1983; BARRet al., 1983; RAMEY, 1985; MEIJERet al., 1987; BODINE et al., 1989; SIEGELet al., 199 1; and cited references). GONZALEZ( 1983) and BARR et al. ( 1983) proposed that groundwater in the Culebra flows southeast from the center of the WIPP site and then westward to Malaga Bend. RAMEY (1985 ) observed that the water chemistries in the Culebra along those hydrologic flow paths require transformation of a saline NaCl brine to a dilute Ca-SO., solution. No plausible geochemical * Paper presented at the “Symposium in Honor of Heinrich D. Holland” held May 8 and 9, 1992, in Reston, Virginia, USA. 2299

process has been identified that would cause this transformation in a hydrologically confined unit. Ramey also defined hydrochemical facies for Culebra waters and calculated saturation indices for common evaporite minerals. The purpose of the present study was to formulate a conceptual model for geochemical evolution of groundwater in the Culebra dolomite during the last 12,000 to 16,000 y. A critically evaluated geochemical database, including analyses of major and minor solutes, mineralogical studies of intact core samples, and isotopic studies of waters, carbonates, and sulfates, was assembled for the Culebra dolomite and adjacent units in the northern Delaware Basin. The large geochemical and hydrological databases available from the WIPP site provide an excellent opportunity to apply several different approaches that have been developed by hydrochemists for the analysis of groundwater chemistry data. These include analyses of spatial distributions of element concentrations and ratios (hydrochemical facies and geochemical signatures), interelement correlations expressed by principal component analysis, and thermodynamic constraints (mineral saturation indices). These relationships were compared to thermodynamic reaction path calculations and models of net geochemical reactions along modem flow paths. Each of these techniques has limitations, requires ad hoc assumptions, and provides multiple possible interpretations of the solute data. One of the objectives of this work was to compare the results of the different geochemical models with geohydrological data from the site, in order to identify an internally consistent set of processes that have strongly influenced the chemical evolution of water in the Culebra dolomite. These same processes

2 300

M. D. Siegel and S. Anderholm

may affect the migration of radionuclides released from the WIPP in the future. GEOLOGICAL

if nuclear waste is

SEITING

Detailed descriptions of the stratigraphy and geological features of the northern portion of Delaware Basin in southeastern New Mexico aregiven by POWERSet al. ( 1978), LAMBERT( 1983), LAPPIN( 1988) and SIEGELand LAMBERT( 199 I ). The Delaware Basin became a distinct structure approximately 250-280 million years ago (Ma) when the marginal Goat Seep Dolomite and overlying Capitan Limestone and the basinal sandstones, shales, and carbonates of the Delaware Moun~in Group were deposited (Fig. t ) . The uppermost of the basinal sandstones and shales (Bell Canyon Formation) grade upward into the alternating sequence of sulfate and carbonate laminations characteristic of the Castile Formation. The Salado Formation overlying the Castile consists dominantly of halite, interrupted at various intervals of meters to tens of meters by beds of anhydrite, polyhalite, mudstone, and local potash mineralization. The Rustler Fo~ation, overlying the Salado, contains the following hve members in ascending stratigraphic sequence: ( I ) a lower member ( unnam~) consisting of sil~tone/mudstone, anhydrite locally altered to gypsum, and containing halite under most of the WIPP site; (2) the relatively permeable Culebra Dolomite Member, containing primarily dolomite with subordinate anhydrite and/or gypsum; (3) the Tamarisk Member, dominantly anhydrite/gypsum, with subordinate fine-mains &sties, confining halite to the east ofthe WIPP site; (4) the Magenta Dolomite Member, another dolomite/sulfate unit: and (5) the Forty-niner Member, similar in lithology to the other non-dolomitic units, but containing halite only in the easternmost part of the study area. The Dewey Lake Red Beds are the uppermost Permian unit, consisting of siltstones and claystones: the Triassic Dockum Group (undivided) rests on the Dewey Lake in the east halfof the site and thickens eastward. Channel and alluvial pond deposits of the middle Pleistocene Gatuiia Formation occur locally throughout the area.

Mineralogical analyses of more than IO0 samples of fracture matrix and coatings from locations along three east-west trending traverses and WIPP shafts have been carried out ( KRUMHANSLet al., I 99 I: &WARDS et al., 199 1a,b,c, 1992; SEWARD&I99 1). The dominant mineral in the Culebra dolomite is a fairly pure dolomite (Ca/Mg = - 1.05; Fe0 = 0.23%) and comprises about 85%, by weight of the bulk rock. TEM and XRD crystallographic studies showed that the dolomite is homogeneous. well ordered, and essentially ideal (Sr;.WARDS,1993). Fractures are observed in all cores and are lined most commonly with clay and gypsum. Clay minerals are present in all samples occurring either as discrete seams, along fracture surfaces or finely dispersed throughout the matrix. They comprise approximately 1to 7% by weight of the bulk matrix and coat 1to 43%’of the fracture surface examined. In the Culebra, the dominant clay minerals, in order of abundance, are corrensite > illite > chlorite :> serpentine. On the regional scale, corrensite, which is an ordered interstrati~ed saponite/chlorite, is the most abundant mineral observed after dolomite. Gypsum occurs as fracture and vug fillings; all fracture surfaces examined by XRD showed the presence of nearly pure CaSO, .2HrO. Calcite is present in significant proportions in the top part of the Culebra samples from WIPP-29 and is interpreted to be of secondary origin ( SEWARDSet al., 199 1a). Microprobe analyses of microcrystalline calcite in core samples from WIPP-I9 show that it is a very low-Mg calcite. Trace amounts of pyrite, magnesite. quartz, and authigenic feldspar have been observed in some cores; occurrences of halite and other evaporite salts are very rare. HYDROLOGIC

SETTING

Studies of the hydrologic system in the Cuiebra dolomite and related units are summarized by LAPPIN( 1988) and SIEGELand LAMEERT ( 199 1) and are briefly reviewed here. Modern regional flow within the Culebra at and near the WIPP site appears to be mainly northto-south. Calculated Darcy velocities range over six orders of magnitude, from IO-” m/s (m3/m2s) east of the WIPP site to as high as 10m6m/s in central Nash Draw (Fig. 2). Estimated values of

FIG. 1. Setting of the WIPP site relative to the northern Delaware Basin. Selected geomorphic features and some boreholes referred to in text are also included. From LAPPIN( 1988) modified from Figure VI-I of LAMBERT( 1983).

2301

Geochemical evolution in a dolomite aquifer

R29E

MOE

I

*Fv

RUE

I

R32E

I I I

.(-8 ‘n’ :

,I’v

M

” *f ,,,

;; I-

4,

.r +

P ..d

,a +

f

-I-

% * s 0

e

LIVINGSTON

.AEC-9

4k 9 Q \‘ &

0

RIDGE SURFAeE

DOE-20 _ _WIPP-14 4 H-5 WIPP-34

H-6

f

I I I I I

I L :

(II

E WIPP-26

. WIPP-32

l

<- I,,.

*H-3

8 l

WIPP-29

*DOE-l

H-14*

i. +

b

t,

-H-4 l P-17* Windmill Well f Re,,ch WeII. c

w

5, Bern We! l

lSouth wetl

FR-10

.H_7

4

$

.

I

llndlen Well f * g,. ,* .Y

.

*H-17

I I

1 I

“0

1

Un9erWetl

*H-12

0 + LIVINGSTON

I I

\”

5. * home-5

I I I I I

Pocket Well Fairview Well*

#ERDA-10

ym *Gnome Sheit

.*Gnome_4 = Gnome-O

+

l

’

’

u

RIDGE hRFACE

%Y ‘4.wi

_. ),SL ....

’ I + I

*H-l1

_ P-15

c. -4. -d

*P-l8

E

H-10.

I

Gnome-l

+-

-I*H-9 l

Two-Mile

Well

l

l

.-.’

Engle Well

Twin W*

r: t

OH-9

,. r

Poker Trep W’,ll

0 R29E

I

R30E

I

R31E

2

4

/

6 km R32E

FIG.2. Map showing locations of wells mentioned in this article. From SIEGELand LAMBERT ( 199 1).

Culebra transmissivity in the WIPP site area and within Nash Draw range over approximately six orders of magnitude, from about 2 X lo+ m2/s a short distance east of the site, to about 1 X low3 m’/s west of the site, on the east side of Nash Draw; in the central portion of the site, estimated values are less than IO-’ m2/s (LAVENUE et al., 1988). The degree of fracturing and fracture filling appears to have some control over transmissivity in the Culebra, especially where the transmissivity is higher than about 10m6m*/s. The amount of possible vertical flow into and out of the Culebra

near site-center has not been directly measured in the field. At several well locations, the relationships between the heads in the Culebra and adjacent units indicate the potential for vertical flow into the Culebra from the Magenta (H-16, H-3, H-2, H-l, H-14, and DOE1) or the Rustler/Salad0 contact zone (H-2 and H-16); however, at other sites (holesH-3, H- 16, and H- 14), the head relationships between the Forty-niner and the Magenta preclude recharge of the Magenta and underlying units by direct downward infiltration from the surface( BEAUHEIM,1987;SIEGELand LAMBERT, 1991). The ex-

2302

M. D. Siegel and S. Anderholm

tremely low hydraulic conductivity of the Tamarisk and the tower unnamed member and relatively high transmissivity of the Culebra suggest that confined flow within the Culebra prevails at and near the WIPP site.

.

,*-

l-‘--.ii._

LAMBERT ( 199 I) suggests that stable isotopic signatures indicate that it is unlikely that the confined Rustler and Capitan groundwaters in the northern Delaware Basin arc receiving modern meteoric recharge under the present climatic conditions. Demonstrably modern

=R

I

i

en17

FIG. 3, ~e~ation~~~ between modern flow inferred fromb numerical modeling of physical hydrology data (solid arrows) and postulated paleofiow (late Pl&%o@xx ) directions fof groundwater in the ~~~.~~m~r of the Rustfei Formation. Lwations of the WIPP !SiteBoundary (solid Tines)and Nash Draw limits (hachunxl Iii)-and hycirw~m&d facies discussed later in Paper are shown for reference, Compositions of waters in each zoneare described in Fig. 5, From SIEGEL and LAMBERT( 199 I).

Ceochemical evolution in a dolomite aquifer

waters in the region, such as vadoseCapitan waters from the Carlsbad Caverns system,groundwater from a water table in alluvium, Ogallala groundwatersunderlying the High Plains to the east, and near-surface accumulations all have meteoric isotopic signatures that are significantly more positive in 6D-6i80 space than confined Rustler and Capitan waters. CHAPMAN( 1988) provides an alternate view and suggests that available data are consistent with modem recharge of the Culebra. Radiocarbon data suggest that confined Rustier and Dewey Lake groundwaters have a significant component at least 12,000 radiocarbon years old. Due to anthropogenic contamination, only four well water samples (three Culebra and one Dewey Lake) that give internally consistent radiocarbon “age” calculations are available ( LAMBERT, 1991). The minimum mean age of these four is 14,000 years; 95% confidence limits are ?5,000 years, based on the interpretive model of EVANS et al. ( 1979) for water-bearing carbonate rocks. LAMBERT ( 199 I ) suggests that zJ4U/z38U relationships indicate that general paleoflow directions in confined Rustler groundwaters involved recharge with an eastward-flowing component (Fig. 3). This is not parallel to modem flow directions inferred From potentiometric heads. He suggests that flow direction has changed back toward the west (the modem flow direction indicated for saturated portions of the Magenta; MERCER,1983) during the past 10,000 to 30,000 years. This is consistent with the results of studies summarized by LAMBERT

( I99 1) which suggestthat the J-atePleistoceneclimate of southeastern New Mexico was probably wetter and more conducive to recharge than the present. During the time of wetter Pleistocene conditions, if the Rustler was being recharged in outcrops in Nash Draw, the base level there was probably being lowered by accelerated erosion and dissolution. Since the onset of drier Holocene climate and cessation of significant recharge, this lowered area has probably governed the local hydraulic gradient, resulting in the general southward drainage observed today. HYDROCHEMICAL

DATA

Sources of Data

Analyses of major and minor solutes from water samples collected by several different agencies and analyzed in different laboratories from 1976 to I987 were examined in this study ( MERCER, 1983; MERCER and ORR, 1979; UHLAND and RANDALL, 1986; UHLAND et al., 1987; RANDALLet al., 1988; CHAPMAN, 1988). In many cases, large discrepancies exist in the concentrations of major and minor solutes reported for the same location by different agencies. Much of the variation can be attributed to differences in sampling methods, to differences in analytical techniques, or to changes in the groundwater chemistry due to contamination. LAMBERTand HARVEY( 1987) and SIEGELet al. ( 1991) describe the different datasets that were examined and the development of an evaluation procedure to determine the criteria for acceptable values. Where available, water samples from the Water Quality Sampling Program ( WQSP; UHLAND and RANDALL, 1986; UHLAND et al., 1987; RANDALLet al., 1988; LYON, 1989) were used for the calculations described in this article. Many of the WQSP samples have been analyzed by up to three different laboratories, and the results are amenable to interlaboratory comparisons. In general, most of the major, minor, and trace solute data obtained at UNC Geotech (formerly Bendix Field Engineering Carp, Grand Junction, CO, USA) are considered reliable and are used as reference values. Table Al lists the values that were used in the trilinear diagram, saturation index calculations and factor analyses; Table A2 lists the values used in the NETPATH calculations discussed below. The compositions in Table A2 were defined after re-

2303

viewing more recent chemical analyses of several samples from each well (ROBINSON, 1992). In some cases they differ slightly from the earlier analyses listed in Table A 1. It is assumed that samples are representative of the given location if samples of the same composition have been collected during several different sampling periods several years apart. Samples obtained using a “serial sampling” method have been shown to be more reliable than those obtained by bailing or swabbing ( LAMBERTand ROBINSON,1984). During a pumping test, the “serial samples” are collected at regular intervals and analysed for concentrations of several solutes. When the chloride and divalent cations (Ca plus Mg) reach steady-state concentrations, the samples are considered representative with respect to concentrations of the major solutes (Cl, S04, Ca, Mg, Na, and K) and certain minor and trace elements (Br, F, I, B, Li, and Sr). At least two wells (H-5, WIPP-25) may have been contaminated with modem organic carbon during the drilling process ( LAMBERT, 199 1); these organics were apparently metabolized by bacteria, leading to high transient alkalinities. Other potential sources of error associated with characterization of the carbonate chemistry of Culebra waters are discussed in a later section. SOLUTE RELATIONSHIPS IN WATERS FROM THE CULEBRA Spatial Distributions of Major Solutes

In simple hydrologic and hydrochemical systems, spatial trends in element concentrations and ratios can be used to suggest flow directions and to delineate potential sites of recharge or vertical leakage. Average concentrations of major and minor elements in Culebra waters are presented as contour plots in SIEGELet al. ( 199 1). Regional variations in the concentrations of Na are contoured in Fig. 4; the actual values for each well are also given. Contour placement is controlled by several factors including the selected contour interval, smoothing factors, and radius of influence of the data points; therefore, there are discrepancies between the placement of contour lines and the concentrations at some well locations. The contours show lower Na concentrations in the southern and southwestern areas (H-7B, H-8B, H-9B, and Engle), and higher concentrations in the north and east. This pattern is also representative of the distributions of K, Mg, Cl, total dissolved solids (TDS), the Mg/Ca molar ratio, and fluid density. Contours of silica concentrations show an inverse pattern. Data for WIPP-27 and WIPP-29 were not used in construction of the contour plots. As discussed in SIEGELet al. ( 1991), the K/Na weight ratios, the Mg/Ca molar ratios, and Na and Sr concentrations are anomalously high at these sites, suggesting that contamination by potash refining operations is important at these locations. Hydrocbemical

Facies

Based on the major solute compositions described in Table A 1, four hydrochemical facies or zones are delineated in Fig. 3 and 5: Zone A. A saline (about 3.0 molal) Nail brine with a Mg/Ca molar ratio between 1.2 and 2.0. This water is found in the eastern third of the site; the zone is roughly

1304

M. D. Siegel and S. Anderholm :oincident with the region of low transmissivity. On :he western side of the zone, halite in the Rustler nas been found only in the unnamed lower member; m the eastern portion of the zone, halite has been observed throughout the Rustler (see Fig. I - I3 in SIEGEL et al..

WIPP WC /

BWIPP-2s

190 O:::: 0

2

4

m’ 6 km

H-6 ?

FIG. 4. Contours of Na concentrations (mg/ L) in the study area. The location of the 4 mile zone is indicated for reference. From SIEGEL and LAMBERT ( 199 1).

1991 ).

Zone B. A dilute CaS04-rich water (ionic strength < 0. I m) occurs in the southern part of the site. Mg/Ca molar ratios are uniformly low (0.4-0.5). This zone is coincident with a high-transmissivity region; halite is not found in the Rustler in this zone. with low to modZone C Waters of variable composition erate ionic strength (0.3- 1.6 m) occur in the western part of the WIPP site and the eastern side of Nash Draw. Mg/Ca molar ratios range from 0.5 to 1.2. This zone is coincident with a region of variable transmissivity. In the eastern part of this zone, halite is present in the lower member of the Rustler; on the western side of the zone, halite is not observed in the formation. The most saline (Na-Cl-rich) water is found in the eastern edge of the zone, close to core locations where halite is observed in the Tamarisk Member. wells WIPP-27 and Zone D A fourth zone, comprising WIPP-29. can be defined based on inferred contamination related to potash-refining operations in the area. Waters from these wells have anomalously high salinities (3-6 m) and K/Na weight ratios (0.22) compared to other wells at the site (O.Ol0.09). At WIPP-29, the composition of the Culebra water has changed over the course of a 7-year monitoring period. The Mg/Ca molar ratio at

(c--v\ /

FIG.

5.

Trilinear diagram for Culebra groundwaters. From

DDE-2

SIEGEL

\

et al. ( 199 I )

2305

Geochemical evolution in a dolomite aquifer

SIEGELet al. ( 199 1) summarize the basic concepts and vocabulary of principal component analysis (PCA) and describe the application of this technique to waters from the Culebra dolomite and related rock units. The principal components identified from Q-mode PCA can be thought of as hypothetical compositional endmembers; principal components obtained from R-mode PCA express the relationships among the compositions in terms of a small number of independent underlying correlations. For example, two cations that are dissolved from the same mineral will be positively correlated whereas elements whose concentrations are inversely related by a solubility product will be negatively correlated. This information can be used to refine the definition of hydrochemical facies and to suggest the nature of chemical reactions that control groundwater composition. Note that in this article, the terms “factor” and “principal component” are used interchangeably. However, the analysis presented below is in the strictest sense a principal component analysis and not a classical factor analysis.

WIPP-29 is anomalously high, ranging from 10 to 30 during the monitoring period. The chemical characteristic3 of each of the facies are described by the Piper (trilinear) diagram shown in Fig. 5. The diagram shows the relative proportions of the ions on an equivalents/ liter basis in the Na-K-Mg-Ca-Cl-S04-COs system. Compositions of Culebra waters from Zone A have nearly identical ionic proportions and plot on the same location in the graph near the Na-Cl comer. They are distinguished from waters in Zone C by their higher Mg/Ca ratio and ionic strengths (Fig. 6). Dilute groundwater from Zone B plot near the CaSO4 comer and water from Zone C cover a wide area in the plot. Waters from Zone D have similar ionic proportions to those of Zone A and are distinguished primarily on the basis of salinity and K/Na ratio, parameters that are not shown in Fig. 5. Principal Component Analysis of Water Chemistry Data The relative proportions of the concentrations of seven major solutes can be easily represented by a trilinear plot as shown in Fig. 5. The method of principal component analysis (PCA) was used to show relationships among thirteen chemical variables (Ca, Mg, K, Na, Cl, SO,, B, Li, SiOz, Br, Sr, alkalinity, and pH) in twenty-one Culebra water samples. This method summarizes the relationships among the variables and shows how these relationships differ for each water composition. The factor analysis programs of SAS (Statistical Analysis System Institute, 1982) were used to extract the principal components. The dataset consisted of the logs of the solute concentrations listed in Table Al. Where multiple samples from the same well are listed, the most recent sample was used; wells H-4C, H-X, H- 11B, and WIPP- 13 were not included in the analysis. All principal components (factors) that accounted tracted.

for at least 1% of the total variance

3

I

I

Results A Q-mode principal component analysis indicated that two principal components (factors) account for nearly all of the variance of the data. In general, if two or three Q-mode principal components account for most of the sample variance, then the population can be considered homogenous enough for R-mode analysis (see HITCHONet al., 197 1; KLOVAN, 1975; DREVER, 1982). The results of the unrotated solutions and more detailed discussion ofthe Q-mode analysis are found in SIEGELet al. ( 199 1). The R-mode correlation matrix and factor loading matrix for the Culebra solute data are given in Tables A3 and A4, respectively. Three principal components (factors) account for approximately 96% of the total variance of data (Table A4 and Fig. 7.) Factor 1, the “salinity factor,” accounts for approximately 67% of the total variance of the dataset. It is

were ex-

I

I

1.50

2.00

o 5 a

0.50

1.00

3.50

4.00

lonlc Strength (mold)

FIG. 6. Relationship between Mg/Ca molar ratio and ionic strength (molal) for Culebra groundwaters. From SIEGEL et al. (1991).

M. D. Siegel and S. Anderholm

NJ

Cl

CJ

Mg

K

Sr

pH HCO,

so,

St02

Bf

B

LI

1.00

1

Na

Cl

Ca

Mg

K

Sr

pH HCO, SO, SiO,

Br

B

Li

CJ

Mg

K

Sr

PH

Br

B

Li

-0.50 Factor 3

-1.00

FIG.

7.

i NJ

cl

HCO,

so,

902

Unrotated R-mode factor loading for Culebra groundwater samples. Values are listed in Table A4. From

SIEGEL et al. ( 199 I ).

characterized by high positive correlations (loadings) for all of the variables except pH, bicarbonate alkalinity, and silica. The factor accounts for most of the variance of Na and chloride, and it is likely that this factor represents addition of solutes by dissolution of evaporite minerals. Factor 2, the “silicate/bicarbonate factor,” accounts for approximately 18% of the variance of the dataset. It is characterized by large loadings for pH, bicarbonate alkalinity,

and silica and B concentrations. This combination of elements may reflect sorption and silicate and carbonate diagenesis as discussed in a later section. Factor 3, the “sulfate factor,” accounts for approximately 11% of the total variance and exhibits high positive loadings for Ca and Sr and a high negative loading for sulfate. This element association probably reflects dissolution/precipitation of gypsum. The unrotated principal components described above were

2307

Geochemieal evolution in a dolomite aquifer compared to those obtained from several alternative PCA rotations (varimax, equamax, and quartimax rotations) to determine if any element associations were common. It was found that the factor-loading matrices from all rotational schemes contained the same key element associations described above for the unrotated components. The loading matrices differed primarily in the relative loadings of other elements and the variance explained by the factors. The compositions of the samples in terms of the factors are described by the factor scores. Relationships between varimax factor scores for the sahnity and sulfate factors described by SIEGELet al. ( 199 1) show that the first two principal components can be used to delineate the same groups of wells identified in Figs. 3, 5, and 6. Thus, the definition of hydrochemical facies based on the proportions of major solutes is supported by correlations among the major and minor elements examined in this study. R-mode factors obtained after partialling-out total dissolved solids The three major factors obtained from the R-mode analysis described above are strongly affected by halite dissolution. Solutes are added directly from the halite or indirectly because the solubilities of sulfate and carbonate phases increase as the ionic strength increases (see next section). A second Rmode principal component analysis was carried out to examine interelement correlations independent of the effects of halite dissolution. In this analysis, the principal components of the partial correlation matrix with respect to total dissolved solids concentration were extracted as follows. First, regression equations for each of the chemical variables as a function of the TDS were obtained. Next, the partial correlation matrix with respect to TDS was obtained by calculating the correlations between the residuals from the regression ~uations. Finally, the eigenvectors of the partial correlation matrix were extracted to give the principal components. A more detailed description of the procedure can be found in SIEGELet al. ( 199 I). Five factors account for 99% of the variance that remained after the TDS was partialled out; the factor loading matrix is found in Table A5. The amount of the variance that is not correlated with TDS is shown for each variable in the last column of Table A5. A significant portion of the variation in the concentrations of Ca, Sr, B, and Siq, the pH and the alkalinity cannot be correlated with the TDS. In Fig. 8, the factor loadings from Table A5 are recalculated in terms of the percent of the total variance (including the variance correlated with the TDS) of each solute that is explained by each factor (note that the variance is equal to 100 X (factor loading)’ X residual variance). Factor 1B is dominated by the negative correlation of Ca and Sr with sulfate and accounts for 67%, 35%, and 24% respectively, of the total variance of these elements. It is suggestive of dissolution /coprecipitation of Sr and Ca in a sulfate phase such as gypsum or anhydrite. Factors 2B, 4B, and 5B contain all or portions of two negatively correlated groups of variables. One group involves the correlation of Mg, K, bicarbonate alkalinity, and silica; the other group contains pH, B, and Li. This pattern of element association is not due to dissolution of evaporite minerals. It

may reflect a combination of processes such as sorption, ion exchange, carbonate diagenesis, and silicate diagenesis, and is discussed in more detail in a later section. Saturation Indices for Culebra Waters Calculations of saturation indices of groundwaters are useful in suggesting the identify of the minerals that may control solute concentrations, in detecting evidence of metastable persistence of mineral phases, and in indicating possible errors in chemical analyses or sampling procedures. Saturation indices (SI) are calculated as SI = log IAP - log K,

(1)

where the IAP is the ion activity product of the species involved in the dissolution reaction, and K is the thermodynamic ~~lib~urn constant for the reaction at the relevant temperature and pressure. At saturation, the saturation index equals zero; positive and negative values indicate supersaturation and undersaturation, respectively. Saturation indices for common evaporite minerals for Culebra brine samples are presented in Table 1. The calculations were carried out using the geochemical reaction and speciation code PHRQPITZ (Version 0.1, unscaled option; PLUMMER et al., 1988). This code uses a specific interaction model to calculate activity coefficients in the system Na-Ca-MgK-H-OH-COz-SO4-Cl. The theory, data, and implementation are based primarily on previous work by Pitzer and coworkers (PITZER, 1973, 1975; PITZER and KIM, 1974; PITZER and MAYORGA, 1973, 1974) and Harvie and coworkers (HARVIE and WEARE, 1980; HARVIE et al., 1984). Estimates of the maximum ionic strength for which the model accurately predicts mineral solubilities range from 6 to 20 molal. Table I and Fig. 9 show that all water samples are saturated f Sl = 0.00 f 0.09) with respect to gypsum f CaSO, * 2HZO). This is expected in tight of the ubiquitous presence of gypsum in the Culebra and the rapid dissolution rates of this mineral under natural conditions (BACK et al., 1983). With the exception of the Culebra samples from WIPP-29, all samples are unde~tumted with respect to anhydrite. All of the waters analyzed from the Culebra are undersaturated with respect to halite. Saturation indices were also calculated for other common evaporite salts; all waters were undersaturated with respect to these minerals (SIEGELet al., 1991). Saturation indices of carbonates Table 1 shows that many of the samples show appreciable supersaturation or undersaturation with respect to calcite and dolomite (Fig. 9). Given the radiocarbon ages for Culebra waters (at least 10,000 years old), the apparent lack of carbonate equilibrium is problematic. The near-zero saturation indices for gypsum suggest that errors in analysis of Ca are not responsible for possible errors in the saturation indices for the carbonates. Other possible sources of error are discussed below. Errors related topfi measurements. PLUMMERet al. ( 1988) show that the magnitude of the error in calculated saturation indices introduced by the use of different activity scales measured pH and computed ion activities in a speciation cal-

M. D. Siegel and S. Anderholm

2308

60

-40

A

Ca

Sr

so,

-40 -60

_-

HCO, SIO,

8

Li

FIG. 8. Amount of total variance of key etements explained by factors 13 to 5B obtained from the ~~i~~o~elation matrix with respect to TDS of Culebra waters. Complete factor-loading matrix is found in Table AS. From SIEGEL et al. (1991). culation can be large for the carbonate system at moderate to high ionic st~ngths (>0.7 molal). In the Culebra samples, however, the degree of disequilibrium is not related to ionic strength (Table I ). Therefore, increase of activity scale error with increasing ionic strength cannot be the sole cause of the apparent lack of saturation. uncertainty in thermodynamic data for ~~nerai phases. Nonstoichiometry and disorder can be a source of uncertainty in the Gibbs free energy of formation of dolomites ( REEDER and WENK, 1979; REEDER, 198 1). However, detailed X-ray diffraction studies, TEM studies, and chemical analyses of dolomite samples from four sites in the study area indicate that the Culebra dolomite is well-ordered, free of exsoiution, and is essentially pure CaMg( COsh (SEWARD&1993 ). This suggests that any error in the calculated saturation indices due to the presence of disordered dolomite in the Culebra is probably not sibilant. Errors related to loss of C?& gas during sampling. Changes in pH due to loss of CO2 gas during water sampling result in apparent supersaturation of carbonate minerals; H+ produc-

tion from oxidation of Fe*+ results in apparent undersaturation ( BASSET, 1982 ) . For solutions at ~uilib~um with eaicite in the pH range of waters in the Culebra (pH -6 to S), the error in the measured pH, ApH (where ApH = pHmeasud - pHequilibri”m),is equal to the saturation index (i.e., ApH = SI; PEARSONet al., 1978 ). SIEGELet al. ( 199 1) demonstrate that some portion but not all of the calculated su~~tumtion is attributable to this effect. The degree of carbonate mineral saturation, independent of the effects due to CO2 outgassing, can be determined by examining the value of the SI expression (PEARSON et al., 1987) 2

%alcite

-

S~dolomiter

(2)

where SI, the saturation index for each mineral, is defined by Eqn. 1. The value of this expression is independent of the unbind in the carbonate ion activity and equals zero when the solution is saturated with respect to both calcite and dolomite. Figure 10 shows that the values of the SI expressions for Culebra waters decrease as a function of ionic strength.

2309

Geochemical evolution in a dolomite aquifer Table

1. Mineral Saturation Indices for Culebra Groundwaters Saturation

Ionic

Temp

2

Date

Well

Indices’*S

Str.l

pH3

&O

2.53

7.10

0.923

CBE1s4

Anh

Cal

Dol

P.OSE-03

-0.12

-0.34

-0.13

Gyp

Hal

Mag

logpC02’

0.03

-1.26

-0.62

-2.60

(“C) DOE-1

4185

22

DDE-2

3185

22

1.19

7.00

0.968

ENGLE

3185

22

0.09

7.40

0.999

H-2A

4186

22

0.27

8.00

0.994

H-3B3

2185

22

1.08

7.40

H-3B3

6184

26

I .08

7.40

-0.18

-0.21

-0.14

0.01

-2.03

-0.77

-2.33

-0.22

0.22

0.38

0.00

-6.02

-0.68

-2.44

-4.338-03

-0.25

0.45

0.78

-0.04

-3.49

-0.52

-3.38

0.971

-2.538-03

-0.20

-0.08

0.12

0.00

-2.11

-0.65

-2.86

0.972

-2.568-02

-0.15

-0.02

0.23

0.04

-2.13

-0.59

-2.83 -3.35

-6.52E-02 I .OOE-03

H-lB

S/8 1

22

0.455

8.00

0.991

- 1.498-02

-0.17

0.36

I .07

0.04

-3.13

-0.13

H-4B

7185

22

0.42

7.70

0.991

4.14E-03

-0.20

0.07

0.46

0.01

-3.17

-0.45

-3.04

H-4C

8/84

21

0.45

7.80

0.991

4.82E-03

-0.20

0.19

0.77

0.01

-3.12

-0.26

-3.11

H-5B

618 1

22

2.993

7.90

0.906

-1.18E-01

-0.10

0.71

2.12

0.04

-1.08

0.57

-3.21

H-5B

8/85

22

2.97

7.40

0.908

8.46E-02

-0.11

-0.02

0.68

0.03

-1.09

-0.15

-2.86

H-K

10181

24

3.002

7.90

0.906

-1.20E-01

-0.08

0.75

2.20

0.05

-1.08

0.61

-3.17

H-6B

S/8 I

22

1.178

7.00

0.969

2.248-04

-0.15

-0.04

0.16

0.05

-2.05

-0.64

-2.16

H-6B

Y/85

22.5

1.13

6.90

0.970

-1.288-02

-0.20

-0.16

-0.07

-0.01

-2.08

-0.75

-2.07

H-7Bl

3186

22

0.089

7.20

0.999

-l.OPE-04

-0.23

0.07

0.01

-0.01

-5.86

-0.90

-2.20

H-8B

1186

22

0.083

7.60

0.999

-1.978-04

-0.23

0.33

0.66

-0.01

-7.45

-0.52

-2.70

-0.22

0.24

0.37

0.00

-6.23

-0.71

-2.43

1 l/85

H-PB

22

0.087

7.40

0.999

3.708-04

H-1183

6185

22.7

2.23

7.20

0.934

-3.60E-02

-0.13

-0.17

0.12

0.03

-1.39

-0.54

-2.63

H-12

8185

24

2.72

7.20

0.916

4.42E-02

-0.12

-0.15

0.34

0.02

-1.18

-0.35

-2.61

P-14

2186

21.5

0.582

6.80

0.988

-9.068-03

-0.15

0.14

0.18

0.05

-3.01

-0.81

-1.81

P-17

3186

21.2

I .67

7.50

0.953

-3.618-02

-0.18

0.10

0.73

0.00

-1.70

-0.22

-2.90

WIPP-25

8180

23

0.26

6.90

0.995

3.218-03

-0.24

0.03

0.04

-0.03

-3.56

-0.84

-1.69

WIPP-26

8/80

22

0.329

6.90

0.993

-5.12E-03

-0.18

-0.04

-0. IO

0.03

-3.36

-0.91

-1.86

WIPP-

8/80

22

0.330

6.90

0.993

-3.7lE-02

-0.23

0.09

0.14

-0.02

-3.32

-0.79

-1.84

I l/85

22

0.371

7.10

0.992

-9.298-03

-0.20

0.10

0.18

0.02

-3.22

-0.77

-2.14

22

2.573

6.30

0.922

-7.628-02

-0.12

-0.38

-0.40

0.03

-I

-0.87

-1.33

266 WIPP-26 WIPP-27

Y/80

WIPP-28

Y/80

22

0.903

6.50

0.976

-2.428-02

-0.26

0.16

0.53

-0.07

-2.27

-0.48

-0.76

WIPP-29

8180

20

4.908

6.10

0.840

-2.06E-01

0.02

-0.72

0.01

0.09

-0.61

-0.21

-0.87

20

4.472

6.10

0.853

-2.748-01

-0.15

-0.45

0.61

-0.07

-0.64

0.21

-1.10

21.5

6.57

5.90

0.773

-2,07E-01

-0.01

-1.25

-0.56

-0.01

-0.20

-0.16

-0.75

8.80

0.986

-2.508-02

I.12

2.37

0.02

-2.73

0.41

in Table

Al

WIPP-

8180

296

12185

.29

WIPP-29 WIPPJO

P/80

‘Unless

otherwise

al., 1988) (atm.),

calculations

0.582

measured

(see footnote

6), the solute data from IJNC

of the ionic strength

and the satoration

2Temperature 3The

21 indicated

(mob@,

the activity

-0.18 Geotech

of water

given

( aH2G),

were used in the PHRQPITZ

the charge-balance

error (see footnote

-4.41 (Plummer

et

4). the log pC02

indices. at the wellhead

during

sampling

test and assumed

in saturation-index

pH values used in the PHRQPITZ calculations; a few differ slightly from those given charge-balance error (CBE in equivalents per kilogram ofwater) was calculated as:

calculations. in Table

AI.

4The CBE 5Mineral &rhese

(eq/kg

H20)

names: samples

= cations

@q/kg

Anh - anhydrite; were analyzed

H20)

- anions (eqikg

Cal = calcite;

by the USGS

H20)

Do1 = dolomite;

Central

Laboratory;

Gyp = gypsum;

Hal = halite;

the data are given

in Mcrccr

This means that when the effects of the CO* outgassing are corrected, many of these waters are either supersaturated with respect to dolomite and/or undersaturated with respect to calcite. The degree of disequilibrium increases with ionic strength; a chemical model that accounts for this trend is described in the next section. REACTION MODELS Summary of Important Solute Relationships

In the previous sections, examination of Culebra solute data revealed the following relationships: ( 1) increases in the

Mag = magnesite. (1983).

concentrations of Na, Cl, Ca, Mg, K, SO., , Br, B, and Li with ionic strength; (2) an increase in the Mg/Ca ratio with the ionic strength; ( 3) an increase in the apparent supersaturation ofdolomite with ionic strength in contrast to the near-perfect saturation of the solutions with respect to gypsum; and (4) a correlation of Mg with SiOr and negative correlation of both these elements with Li and B (the “silicate/bicarbonate” principal component). Several types of processes are likely to affect the water chemistry in the Culebra: ( 1) precipitation and dissolution of carbonates, sulfates, halite, and other evaporite salts; (2) ion exchange and silicate diageneais involving clays that occur

2310

M. D. Siegel and S. Anderholm

-0.5 -0.5

I -0.4

1

I

-0.3

-0.2

+

1 -0.1

0.0

0.1

f

I

I

0.2

0.3

0.4

0.5

Gypsum Selurrtlon Index

FIG. 9. Relationship between gyusum and dolomite saturation indices in Culebra groundwater samples. From et al. (1991).

as fracture linings and matrix inclusions; and (3) physical processes such as diffusion and leakage of solutes into the Culebra from adjacent strata. Two geochemical reaction models were used to evaluate the potential importance of these processes for the Culebra waters. The PHRQPITZ code was used to simulate the effect of mineral precipitation/dissolution on groundwater composition in model systems dominated by dissolution of evaporite salts. The NETPATH code ( PLUMMER et al., 199 1) was used to calculate the net geochemical reactions (ion exchange. precipitation, and dissolution) consistent with changes in water composition along modern flow paths. These two computer codes provide complementary constraints on proposed reaction path models. NETPATH can be used to identify all possible combinations of reactions that are con-

0.2

I

$ + +

z E 4S

sistent with differences between two water compositions; PHRQPITZ applies thermodynamic constraints to test the feasibility of a specific reaction path model. Salt Dissolution: An Irreversible Process in Partial Equilibrium Systems Saturation indices for groundwaters are often used to identify minerals that may control solute concentrations and to determine if the water is capable of dissolving additional rock. However, even if the saturation index of a mineral is zero, the mineral may be continually responding to an irreversible process causing dissolution or precipitation. In this case, the system is said to be at partial equilibrium, and the mineral/water equilibria are shifting along a reaction path

I

I

+ P-14 Run8 1-7

1 ++t -w

0.0

SIEGEI

-0.2

5 (

-0.4

z ” 3

-0.6

+ WIPP-27

z -0.8

-1.0 0

2

1

lonlc

3

4

Strength (molel)

FIG. 10. Relationship between SI expression (2&e,, - &,rollllte) and ionic strengths for C&bra groundwaters. The saturation indices and ionic strengths are given in Table 1. The Sl expression equals zero for waters in equilibrium with calcite and dolomite. Solutions compositions calculated by reaction-path simulations are indicated by solid lines. From SIEGEL et al. ( 199 I ).

2311

Geochemical evolution in a dolomite aquifer

where log Ksdgypjequals the equilibrium solubility product for gypsum at the temperature and pressure of interest calculated from thermodynamic data; ycaz+, yso:-, the activity coefficients for Ca2+ and SO:-, respectively, calculated using the PHRQPITZ code; &o, activity of water, calculated using the PHRQPITZ code; QerP, mcamso, at equilibrium with gypsum for a given solution composition, and mea and mso, are molal concentrations of Ca2+ and SO:-, respectively. Relationships between ionic strength and the apparent solubility products of dolomite and calcite in Culebra waters are similar to that shown in Fig, 11 (SIEGELet al., 1991). Thus, even if the saturation indices of these minerals are zero, dissolution ofgypsum, dolomite, and calcite could occur and lead to increases in the concentrations of Mg, Ca, and SOi- ) if the ionic strength increases. The ability of this process to account for the observed solute relationships in the Culebra is examined by the reaction path modeling described in the next section.

of the following assumptions were considered in the various simulations: ( 1) closed- vs open-system with respect to CO? gas at two partial pressures; (2) identity of accessory evaporite salts used as sources of Mg, K, and SO, ; and ( 3) degree of saturation of the waters with respect to dolomite. The simulations are described in Table 2. All of the simulations contained two parts. In the first part, the compositions of recharge water (referred to as the “starting solution”) in equilibrium with calcite, gypsum, dolomite, and various amounts of CO* were calculated. In Runs 1, 2, and 5- 11, the total amount of CO2 was fixed and the system was closed to CO2 gas exchange. In the second part of the simulations, the addition of solutes (Na, Cl, Mg, Sod, and K) to the starting solution by dissolution of halite and other evaporite salts was simulated. Potential sources of Mg and SO1 for the Culebra groundwaters are polyhalite ( K2Ca2Mg( SO4)2 - 2H20), camallite (KMgClj6H20), and magnesite ( MgCOs ) . All of these minerals have been observed in at least trace amounts in the Rustler. Throughout the simulations, these salts and halite were added to the solution in constant proportions. The ratios of camallite and polyhalite to halite were set by the observed molar KC1 ratio of Culebra waters (about 0.007). The proportion of MgCO3 added (magnesite/halite = 0.02 - 0.025) was set by the amount of Mg added by other salts and the observed Mg concentrations. The progress of the dissolution reaction was specified by adding the halite in successive increments of 0.3 or 0.6 moles until a total of 6 moles of NaCl had been added. The net amounts of solutes (Mg, Ca, CO,, and S04) added or removed from the solution were the amounts needed to maintain saturation levels specified for calcite, gypsum, and dolomite, and the pco,.

Partial Equilibrium Reaction Path Model of Culebra Waters

Results

driven by the irreversible process (PLUMMER, 1984). It is likely that dissolution of halite in the Rustler Formation acts as an irreversible driving force for changes in water composition. Halite and other evaporite salts are strongly undersaturated in groundwaters in the Culebra dolomite (see Table 1). The dissolution of evaporite salts will add solutes directly and also indirectly by increasing the solubilities of gypsum, calcite, and dolomite. Changes in the apparent solubility products of gypsum as a function of ionic strength for Culebra waters are shown in Fig. 11. The apparent solubility product for gypsum in each groundwater was calculated as log

Qwp= 1% &gyp) - bit (w+Yso:-dI*d,

(3)

Method

The details of the calculations are presented in SIEGEL et al. ( 199 1); the most significant results are summarized here. The changes in pH, alkalinity, and pcoz as a function of reaction progress are compared to measured values in the Cu-

The PHRQPITZ code ( Version 0.1, unscaled option) was used to calculate the compositions of waters that would be produced by several hypothetical reaction paths. The effect

-2.40

0

&

-2.60 E

-I

r) .

-3.20

$

’

-I

Zone A

-I---

.

t

kA t

K a”

.

.

Zone

“-

“.“‘------I

A

,I- i Zone _ C

-3.60 0

I

1

I

I

I

I

1

2

3

4

6

6

7

Ionic Strenglh (molrl) FIG.

From

11.Relationship between apparent solubility product for gypsum and ionic strength in Culebra water samples. et al. ( 199 1).

SIEGEL

M. Il.

2111

Siegel and S. Anderholm

‘fable 2 Summary of Parameters for Reaction-Path Calculations log pCOz’ (atms)

RUIl

I

-2

(closed)

Minerals at Fixed Saturation Index (Sl)2,J

Debra in a series of plots in SIE:GliL. et al. ( I99 1). There is no apparent relationship between the pH, alkalinity. and pcoZ measured in the Culebra and the theoretical trends predicted by simple partial equilib~um models. This is consistent with the possibility that the pH values of many of the Culebra samples were affected by loss of CO2 gas during sampling, oxidation of ferrous iron from the well casing, or bacterial production of bicarbonate associated with contamination of the wells by organic compounds ( LAMBERT, 199 1). Run I assumes a simple partial equilibrium mode1 in which pure NaCl is added to the solution while the groundwaters are maintained at equilibrium with respect to calcite, dolomite, and gypsum. As the ionic strength increases. the solubilities of calcite, dotomite. and gypsum change and dedolomitization (dissolution of dolomite and gypsum and precipitation of calcite) occurs, leading to changes in the concentrations of Ca. Mg, and SOa. Figures 12 and 13 show that Run 1 does not reproduce the compositions of Cutebra waters. This simple model provides a good fit to many of the observed calcium concentrations; however, concentrations of SO4 (and Mg) are considerably below the predicted values. It was not possible to reproduce the relationships among Ca, Mg. Cl. and SO4 observed in Culebra waters by dissolution of polyhalite, carnallite. magnesite, and halite if a constant saturation level was maintained with respect to calcite. gypsum, and dolomite (cf. Runs I-9 ). Figure I2 shows that any constant dolomite saturation index leads to a MgiCa ratio

Molar Ratio oi Accessoly Mineral~/~jalit~

c(O).e(O),d(O)

2

-3 I, (closed)

c(O),g(O).d(O)

3

-2 (open)

c(O),g(Oi,d(O)

4

-3,s (open)

c(O).g(O).d(O)

5

-2 (closed)

c(O),g(O),d(O)

poly (0.007)

6

-2 (closed)

c(O),g(O).d(O)

mag (0.02) cam (0 014)

7

-2 (clasedl

401, g(O),d(O)

Icon (0.007)

8

-2 (closed)

c(O),g(O).d(O 5)

car” (0.007) poly (0 007)

9

-2 (closed)

c(O). g(O), d(l.0)

10

-7. (closed)

c(O), g(O)

II

-2 (closed)

m, s(O)

none

leon (0 007) MgC12 (0.025) leon (0.007) mag (0.025)

‘Closed systems were equilibrated with gypsum, calite and dolomite at indicated pC0, beforeadditionof accessory minerals and then closed to iiuther atmospheric exchange. Open systems we% open to exchange of CC+ with ~mosphe~ during entire simulation. ‘Minerals:

c = calcite; g = gypsum; d = dolomite

)Saturation index (SI) for mineral = log (IAP) = log (Ksp) 4poly = polyhalite: mag = magnesite. cam = camallite; leon = leonite.

5

----____ _,_&_ ,/j /

.i”;-,--__-____~~~_~=~~__,

. WtPP-27

1.5

2.0

2.5

lonlc Strenglh fmolrl) FIG _ 12. Change in Mg/Ca molar ratio as function of reaction progress,types of added s&s, and dobmite saturation index (SI) for simulated evolution of Culebra groundwater. Mg/Ca molar ratios of Culebra groundwaters are plotted for comparison. From SIEGEL et al. ( 199 I ).

2313

Geochemical evolution in a dolomite aquifer

I 0

/”

,__,.__. Iy.

/”

+ WIPP-22 mmoo)

....“.““.‘.-‘..“‘.-...

.._. _‘,‘R”nO: b ,;,,

o

_._.”

1

I

I

I

I

I

I

1

2

3

4

5

6

Moles Chloride/kg Hz0

FIG. 13.Change in sulfate concentration as function of reaction progress, types of added salts, and dolomite saturation index (SI) for simulated evolution of Culebra groundwater compositions. Reaction progress is indicated by moles of chloride added. Sulfate concentrations of Culebra groundwaters are plotted for comparison. From SIEGEL et al. ( 1991).

that decreases slightly as the reaction progresses (Runs 1 and 8). This differs from the increase in the Mg/Ca ratio as a

function of ionic strength that is observed in the Culebra waters (Fig. 6). Adequate agreement was only achieved when saturation with dolomite was not maintained (Runs 10 and 1 I ). In these runs, the starting solution was produced by equilibrating distilled water with calcite, dolomite, and gypsum at an initial pcoz of lo-* atmosphere. The system was then closed to COz gas, and addition of Na, Cl, K, Mg, and SO4 by congruent dissolution of halite and incongruent dissolution of polyhalite and camallite was simulated. Camallite dissolves incongruently to form KC1 and a solution containing Mg*+ and Cl-; in Run 10, this was simulated by adding MgC12 to the reaction solution. Incongruent dissolution of polyhalite was simulated by addition of leonite (KzMg(Sod)* - 4H20) and gypsum. The K/Cl ratio of the Culebra waters was used to estimate a reasonable leonite/halite ratio. The MgCl,/halite ratio was then adjusted to fit the observed Mg/CI ratio of the Culebra waters. Equilibrium between the waters and gypsum and calcite was maintained; as the reaction proceeds both calcite and gypsum dissolve. The saturation index of dolomite was calculated as the reaction progressed, but dolomite was not allowed to dissolve or precipitate. In Run 11, magnesite (MgCO3) was added instead of MgClz to adjust the Mg/Cl ratio. Figures 12 and 13 show that Run 10 reproduces the trends of the Mg/Ca ratio and SO, concentrations in this suite of Culebra groundwaters. The relationships Cl, Ca, and Mg in Culebra waters also agree well with the results of Run 10 (see Figs. 2-5 1 and 2-52 in SIEGEL et al., 199 1). The calculated dolomite saturation indices for Run 10 were used to calculate values of the SI expression (2SI,l,it, - SIdolomie)described Previously; the good agreement between the theoretical and observed values are seen in Fig. 10. These calculations show that addition of solutes (Mg, Sod, K) to the Culebra from evaporite minerals such as polyhalite is consistent with the observed groundwater compositions if dolomite does not precipitate from supersaturated solutions.

The reaction-path model described above corresponds to a hypothetical flow model in which recharge water enters the Culebra and dissolves dolomite, gypsum, and calcite. As the water migrates within the aquifer, it obtains solutes from fluids that have dissolved halite and other evaporite salts from adjacent strata. This model is not constrained by the available hydrologic data from the study area, makes ad hoc assumptions, and may not be the only model that fits the observed solute relationships. For example, there is no reason to assume that the ratio of halite to accessory minerals is constant over the entire flow path. Although Run 10 provides the best average fit to the observed Culebra compositions over the entire range of ionic strengths, some of the Culebra waters more closely match the compositions of other runs. Thus, it is likely that each of the water samples in the Culebra has a unique history that does not conform to the simple model described above. To address these issues, the NETPATH code was used to identify other processes that may affect the compositions of Culebra water. Calculation of Net Geochemical Reactions Along Probable Flow Paths Introduction

Figure 14 shows modem flow directions within the WIPP Site ( LAVENUEet al., 1990) from the vicinity of the Air Intake Shaft to the H- 17 well. The salinity increases along similar paths; the Na concentration increases from 14,000 mg/L at the AIS to 54,000 mg/L at H-17. The H-3 and H-l 1 wells are located between the AIS and H-17 and yield waters of intermediate salinity. Estimates of the travel times between wells, based on a steady-state regional flow model for the Culebra ( LAVENUE et al., 1988) range from 2063 years (H11 to H-l 7) to 4225 years (H-3 to H-l 1 ), if the effective porosity of the Culebra is assumed to be approximately 0.08 (SIEGEL, 1992).

The computer program NETPATH (PLUMMER et al., 1991) was used to calculate the mass transfers of phases

M. D. Siegel and S. Anderholm

2314 l

DOE-2

Bounds1

WtPP-She

H-6 H-5

‘WIPP.13 N .WIPP-12

4

element in solution, and J is the number of elements. All possible net geochemical reaction models consistent with the data and candidate phases arc calculated by NETPATH: therefore, individual models must be evaluated against geological data or tested against thermodynamic constraints calculated with PHRQPITZ or similar codes.

lWIPP-16 .WIPP-19

Description

;wlPP.22 (14000) H-16 A,,,Meke

lW’PP-21

q

Shefi(*W

\

.ER,,A.S l

H-15

(81000) Tnvel

Peth

steedy-state Callbnted Field

‘-15 l

H-4 Cabin

Baby-1

i H-17 (54000) 0 t 0

2ml I ’ 3km

1 I 1

2

FIG. 14. Locations of wells considered in NETPATH calculations. Predictedfastflowpaths at the WIPP Site for steady-stateand transient calibrated fields (after LAVENUE et al., 1990;NOVAK, 1992) and Na concentrations (mg/ L) for selected wells are also shown.

(minerals and gases) consistent with the observed changes in the groundwater chemistry between pairs of wells along these flow paths. Given the total concentrations of each element in the initial and final waters of each pair and a list of candidate phases, the NETPATH code writes mass and charged-balanced reactions of the form initial Solution + Reactant Phases = Final Solution + Product Phases,

(4)

where reactant phases and product phases enter and leave the aqueous phase, respectively. The candidate phases can include minerals, gases, or ion exchange reactions. The mixing of two endmember waters in unknown proportions can also be included in the mass transfer to form the final water. NETPATH simultaneously solves a set of mass balance equations for each element of the following form 5 apbp.k = mT.k(tinal)p=,

mT,k

(inittal)

=

Am,k

(5)

(element k = 1 to J.), where P is the number of total reactant and product phases in the net reaction, (Yeis the calculated mass transfer of the pth phase, bp,k is the stoichiometric coefficient of the kth element in the pth phase, mT,kis the total molality of the kth

01’wuter samples

and phases

The mass transfer reactions for the candidate phases considered in the calculations are shown in Table 3. Reactions are written to represent positive mass transfers (i.e., dissolution reactions are positive mass transfers; a positive mass transfer for Ca/ Na ion exchange involves uptake of one Ca ion from solution by the clay and release of two Na ions from a clay to the solution). A hypothetical phase “BRSINK” was included in the calculations to account for mass transfer of bromine. The reference compositions of the waters considered in the NETPATH calculations are described in Table A2. NETPATH calculations are very sensitive to charge imbalances in the water compositions; therefore, the reference compositions in Table A2 were adjusted slightly to eliminate charge imbalances as described in Appendix B. The PHRQPITZ code (Version 0.2, MacInnes Scale) was used to calculate the final compositions in concentration units of mmol/kg Hz0 as shown in Table 4. In the calculations, it is assumed that the composition of the water sampled at a given well (final solution) is produced by net geochemical reactions or mass transfers acting on water (initial solution) that originated near the well lying immediately up gradient in the flow path. For example, it was assumed that water similar in composition to that sampled at H- 18 has flowed to the Air Intake Shaft (AIS), changing composition along the flow path; water similar to that sampled at the AIS has flowed to H-3, etc.

Table 3

Chemical Formulae and Mass Transfer Reactions for Minerals Used in Mass Transfer Modeling’

Halite

NaCl --f Na+ + Cl.

Gyplanh

CaSO‘$ + ca2’

+ so :-

Calcite

C&O3

+ CO ;--

Dolomite

CaMg(CO3)2

+ Ca2’

-_) Ca2’- + Mg2+ + 2C0

;-

Ca@k Ex.

Na2-clay + Ca2+ --t 2Na+ + Ca-clay

Polyhalite

Ca2K2Mg(S04)4 --t 2Ca2+ + 2K+ + Mg2+ + 4S0 :-

Sylvite

KCl+K+tCI-

Camallite

KM&13 --t K+ + Mg2+ + 3CI -

Magnesite

M&O3 + Mg2+ + CO ;

MgMa

Na2-clay

co2

Ex. Gas

co2

+ Mg2+ + 2Na+ + Mg-clay

gas + COZlaa,

K/Na Ex.

Na-clay + K+ + K-clay + Na+

BRSINK

Br

yReactions are written as anhydrous salts because NETF’ATH a constant mass of 1 kg of water.

~~wmes

2315

Geochemical evolution in a dolomite aquifer

Table 4. Water Compositions’ and NETPATH Geochemical Reactions2 Along Modern Flow Paths in the Culebra Dolomite Well

AIS

1 H-3

1 H-11 1

-28.72 13.02 11.85 144.86

0.03 18.71 2.86 22.12 1048.97

11.66 816.57

3.31 256.20

1070.89

0.03

0.28

-0.90

1

l.TizGJ

H-17

1

A~17 - 13-11

-0.04 10.00 3.23 22.62 672.18

9.70

c

A~11 - ~-3

0.36

1

A~17-11-11

]

A

PhaseZ Halite Gypsum/anh Calcite

AHII - ~-3

0.73 49.89 36.16 33.66 766.10

0.36 Mass Transfers’

1 &-MS

256.20 -35.33 -0.43

I

1070.89 -0.68 -12.46 6.31 4.85

1 Water compositions and mass transfers in mmol/kg H20. Compositions were obtained by speciation calculations using PHRQPITZ from compositions in Table A2 followed by a correction for charge balance as described in Appendix B. 2

Reactions defined in Table 3. They are written as anhydrous salts because NETPATH assumes a constant mass of 1 kg water. Positive mass transfers correspond to dissolution of a solid; a negative ion exchange, e.g. Ca/Na, means uptake of 2 moles of Na+ from solution and release of 1 mole of Ca *+ from a solid. Negative CO, mass transfers correspond to outgassing; values are uncertain due to uncertainty in pH measurementsdiscussedearlier in text.

Results

Between ten and fifty models were obtained for each set of wells examined in this study; the large number of candidate phases (Table 3) yields many combinations of phases (mass transfer models) that satisfy the mass balances constraints. Table 4 presents selected net geochemical reaction models for several pairs of wells. A more complete listing of the calculated models is found in SIEGEL ( 1992). As discussed above, the models calculated by NETPATH are constrained by neither thermodynamic principles nor geological data. Several criteria were used to screen the NETPATH models. Models were rejected that ( 1) required precipitation of halite, sylvite, polyhalite, camallite, or magnesite because calculations with PHRQPITZ indicated that all of the waters considered in the study were undersaturated with these minerals (see Table 1); (2) required the precipitation of dolomite, based on the slow kinetics of dolomite precipitation indicated by the reaction path modeling described above; (3) did not include dissolution of halite or precipitation/dissolution of gypsum because of the observed widespread occurrence of these minerals in the Rustler Formation; (4) required dissolution of large amounts of camallite and polyhalite because mineralogical data ( SEWARDS et al., 199 la; MERCER and SNYDER, 1990a,

1990b) suggest that these minerals are very scarce in the Rustler in this part of the WIPP Site; (5 ) required the dissolution of large amounts of magnesite because of the absence of magnesite in the Culebra and its scarcity in other Rustler units; (6) required precipitation and dissolution of large amounts of calcite since available mineralogical data indicate that calcite is rare in Culebra fractures. (In some cases, however, such models were chosen to avoid models requiring Na release by ion exchange as described below); and (7) required release of large amounts of Na from clays in exchange for Ca or Mg from Na-rich solutions because detailed mineralogical studies of Culebra dolomite clays ( SEWARDS et al., 1992) show that the clays are predominantly Mg-rich corrensite rather than Na-rich clays. The remaining models could be grouped into subsets based on the direction and relative amounts of mass transfer of important reactions. “Important” reactions included ion exchange, carbonate phase precipitation/dissolution, and gypsum precipitation/dissolution. Uncertainties in total carbon calculated with PHRQPITZ result from uncertainties in pH measurements and choice of activity scale as discussed previously; therefore, no attempt was made to interpret the mass transfers for CO2 gas. The mass transfer of phases containing the same element are interdependent; therefore, many of the

2316

M. D. Siegeland S. Anderholm

models within a subset differ only by slight differences in the phases included or the amounts of mass transfer. Each model presented in Table 4 is representative of related models within a subset. H-18 to the Air Intake Shaft. If water sampled in the Air Intake Shaft came from and once had the same composition as water sampled in the H- 18 well, reactions resulting in large increases in dissolved Na and chloride concentration have occurred during fluid transport. The increase in Na is larger than the increase in chloride. Sulfur also increases from H18 to the Air Intake Shaft; Ca decreases between the initial and final well. Thirty-nine mass transfer models are consistent with the compositions shown in Table A2; however, all of the models were rejected because they required large mass transfers of calcite, magnesite, or polyhalite, dissolution of dolomite, or release of Na and uptake of Ca by ion exchange. The failure to derive a net reaction path model consistent with conversion of the water from H-18 to water from the Air Intake Shaft may indicate that the two locations do not lie on the same flow path. Air Intake Shaft to H-3. Sodium and chloride concentrations increase between the Air Intake Shaft and the H-3 well (Table 4). The increase in chloride is almost twice the increase in Na. There is a decrease in the S and an increase in both Ca and Mg concentrations. In the fourteen possible models, the increase in Na and chloride was modeled as dissolution of halite. It was assumed that precipitation of gypsum is responsible for the decrease in S concentration. The precipitation of gypsum requires Ca concentration to decrease; however, Ca actually increases between the initial and final wells. The increase in Ca concentration has been modeled by release of Ca and uptake of Na by clays (negative Ca/Na exchange ) . The model shown in the table requires a relatively small amount of calcite precipitation and relatively large release of Mg and uptake of Na by clays (negative Mg/Na exchange). Polyhalite is the source of K, another model involving release of K and uptake of Na (negative K/Na exchange) is also possible. H-3 to H-l 1. Sodium, chloride, S, Mg, Ca, and K concentrations and the Cl/Na ratio increase as water flows from well H-3 to H- 11. Thirty-seven possible models are consistent with the chemical changes. Sulfate concentration can increase due to dissolution of gypsum or polyhalite. If the source of S is modeled by the dissolution of polyhalite, gypsum precipitation is required and either camallite dissolution or ion exchange is necessary. The models described in the table require ion exchange. In model A, dolomite dissolves, a large amount of calcite precipitates and additional Mg is released to solution by ion exchange for Na (negative Mg/Na exchange). Model B is considered less realistic; it requires an unlikely release of Na and uptake of Ca by clays. H-I I to H-l 7. Sodium and chloride concentrations and the Cl / Na ratio increase significantly between H- 11 and H17; concentrations of Mg, K, and S also increase. The main differences among the forty-one possible models are in the amount of polyhalite that dissolves, the amount of gypsum that dissolves or precipitates, the amount and type of ion exchange, the source of K and the identity and amount of carbonate phase which dissolves or precipitates. The models

listed in Table 4 require minor amounts of calcite dissolution, polyhalite dissolution, gypsum precipitation, and uptake of Na in exchange for Mg. Model A requires Ca/Na exchange whereas Model B requires dissolution of carnallite. Possible Changes in Groundwater Composition Due to Mixing of Waters. In the mass transfer calculations discussed in the previous sections, the sources of solutes were limited to reactions involving solid phases (i.e., precipitation or ion exchange). The NETPATH code was also used to compute the mixing ratios and accompanying mass transfer reactions involving both the mixing of selected endmember waters and the reactions described in Table 3 to produce other waters in the Culebra. The mixing of different waters in the Culebra dolomite and the mixing of waters from the Rustler-Salad0 contact zone with water from the Culebra are discussed in detail in SIEGEL( 1992). Many of the calculated models can be rejected because they require large mass transfers of calcite, polyhalite, or magnesite. The remaining models all require mass transfers of phases previously considered in the calculations where mixing was not modeled. This implies that mixing of water masses within the Culebra alone could not produce the observed pattern of water compositions.

DISCUSSION:

COMPARISON OF MODEL AND GEOLOGIC DATA

RESULTS

Several different modeling approaches have been used in this paper to examine relationships among groundwater chemical variables. In the following sections, the consistency among the results from the different geochemical model approaches is described and the proposed processes are evaluated in light of geohydrological data from the site. Dissolution of Evaporite Minerals Dissolution of halite is the dominant reaction in all of the geochemical models described above. However, mineralogical studies ( SEWARDSet al., 1991a,b,c) indicate that halite is relatively rare in the Culebra dolomite within the WIPP site. It is abundant in the Tamarisk and unnamed lower member of the Rustler where it occurs as matrix cement and discrete layers in the bedded mudstone /anhydrite layers. Dissolution of carnallite, polyhalite, magnesite, and/or sylvite is suggested by some of the geochemical models. These minerals have not been observed in the Culebra dolomite, suggesting that the sources of Mg, K, and several other solutes for Culebra waters lie in the layers adjacent to the Culebra. Small amounts of magnesite in halite and mud/&stone matrix or embedded in halite crystals have been observed with X-ray diffraction and scanning electron microscopy in the WIPP-19 core in the unnamed lower member of the Rustler Formation (SEWARDSet al., 199 lc). The presence of the polyhahte has been inferred from gamma-ray logs in core from the lower unnamed member at H- 16, P- 10, and P-l 8 ( MERCER and SNYDER,1990~; SNYDER,1985) and in the Tamarisk, in core from the H-17 well ( MERCER and SNYDER, 1990d). Polyhalite inclusions are common in halite in the Salado Formation (STEIN and KRUMHANSL, 1986). Neither carnallite nor sylvite have been observed in Rustler samples from the Rustler Formation in the vicinity of the WIPP Site.

Geochemical evolution in a dolomite aquifer The mechanism r~~nsible for transport of the solutes from the adjacent layers into the Culebra cannot be determined from the available isotopic, solute, and hydrologic data. Available hydrologic data do not preclude the possibility of relatively slow migration of solutes across the interface between the Culebra and adjacent layers over long f 10,000 year) periods. D/H relationships in gypsiferous horizons adjacent to the Culebra are consistent with the escape of some water from the Culebra into the adjacent beds (LAMBERT, 1991). Dissolution of halite cement or layers in the Tamarisk and lower unnamed member may have occurred simul~eously with gypsification of anhydrite in these beds. Leakage of the water back into the Culebra may have occurred in response to a change in head relationships between the Culebra and adjacent layers. Alternatively, solutes may have diffised from the pore fluids in the Tamarisk or unnamed members into the Culebra. The fluids that dissolve halite cement or layers within the Tamarisk or unnamed members would probably be more saline than water in the Culebra dolomite; the difference in salinity could supply the chemical gradient required for solute di~sion. Precipitation/Dissolution

of Sulfates and Carbonates

Mineralogical studies and geochemical model simulations indicate that gypsum p~ipi~tion/di~olution is an important reaction in the Culebra dolomite. All of the waters are saturated with respect to gypsum, and this mineral will precipitate or dissolve if the Ca or sulfate concentrations change along a flow path. This is consistent with the results of the principal component analysis and is confirmed by the reaction path model. The slight change in gypsum solubility with ionic strength leads to a moderate positive correlation between Na, Cl, and Ca seen in the “salinity” factor obtained in the PCA. This factor accounts for about 35% of the variance of Ca. The equilib~um with gypsum, however, produces a strong negative correlation between Ca and SO, that is expressed in the “sulfate” factor and accounts for more than 50% of the variance in Ca. The negative correlation between pH and bicarbonate alkalinity shown in the “sili~te/bi~nate” factor of the PCA is consistent with dissolution of calcite observed in Run 10: CaC03 + H+ ---t Ca*+ + HCO;.

(6)

The mass action expression for these reactions can be written aS: pH = -log @,coJ - log eaz+ -t log K6,

(7)

where x6 is the equilibrium constant for the reaction described by Eqn. 6. Although dissolution of dofomite may occur locally in low ionic strength waters, the NETPATH calculations indicate that a significant amount of precipitation or dissolution involving carbonate species is not needed to account for the changes in solute concentrations along modem flow paths in the Culebra. Although the solubility of dolomite increases with ionic strength, the supply of Mg from other sources (ion exchange or evaporite salts) prevents dissolution of dolomite; in fact, the PHRQPITZ calculations indicate that the more saline Culebra waters are supersaturated with respect to do-

2317

lomite. Minemlo~~aI data support this hypothesis: the stable isotopic composition of Culebra waters and carbonate minerals indicate no precipitation of dolomite in equilibrium with the groundwater and secondary calcite is rare (SIEGEL and LAMBERT, 1991). Ion Exchange, Sorption, and Silica Diagenesis In the mass transfer calculations for flow between H-3 and H- 1I, and H- 11 and H- 17, Mg/Na exchange (release of Mg and uptake of Na) is required in several ofthe possible models. Such exchange may correspond to the alteration of Mg-rich corrensite to sodium smectite, leading to release of Mg from the clay and uptake of Na from the groundwater. Corrensite, a mixed chlorite smectite, has been observed in Culebra fmctures and in a layer of the lower unnamed member in contact with the Cuiebra ( SEWARDSet al., 1991a,b, 1992; SEWARD& 1993). The clays could be either within the fractures within the Culebra or in the mudstone layers that stratigraphically bound the Culebra. Factor 2B obtained from the partial R-mode principal component analysis of Culebra samples (Table A5 ) shows a moderate Li-B association negatively correlated with a MgSi association. The following amounts of the total variance of these elements are explained by the correlation: SiO,, 53%; B, 35%; Li, 15%; Mg, 6%. The relationships among these elements may be partially controlled by reactions involving silicate minerals. The correlation of Mg and Si in Factor 2B could result from dissolution of a Mg-Si-rich phase such as the trioctahedral Mg-rich layer in corrensite: MgSiX,,,,

+ 2H20 + 2H+ --t Mg*+ + Si(OH)d.

(8)

Evidence that an amorphous Mg-rich layer (MgSiX,,,,) in authigenic mixed-layer clays is more reactive than the AlFe-rich smectite layer is found in studies of Great Salt Lake sediments by SPENCERet al. ( 1985 ) and sediments in alkaline Lake Albert (JONES and WEIR, 1983). High concentrations of clays in fractures could locally lead to elevated concentrations of Mg and Si due to clay dissolution or low concentrations of Li and B due to sorption of these elements by clays. Figure 15 shows that the Li/Cl ratio decreases as the Li concentration increases over most of the ionic strength range. The depletion of Li relative to Cl in these waters could be due to increasingly lower Li contents of the salt that supplies the Cl, or removal by clays of Li supplied by the salt dissolution. The second mechanism is consistent with the wetldocumented uptake of Li by clays in saline Li-rich waters (HOLSER, 1979; COLLINS,1975; G. SMITH, 1976; C. SMITH, 1976) and from seawater ( SEYFRIEDet al., 1984). About 43% of the variance of B is not correlated with the total dissolved solids concentration (cf. Table A5 ). Like Li, B can be either released or sorbed by clays in saline waters. There is abundant field and laboratory evidence that illite and montmorillonite effectively scavenge B from saline solutions such as seawater ( SEYFRIEDet al., 1984, and cited references; HARRIS& 1969; GOLDBERGand GLAUBIG, 1986; SONNENFELD,1984; HARDER, 196 I ) . The distribution of aqueous silica at the WIPP site may be controlled by saturation with respect to a silica or silicate phase, the rate of dissolution from a silica-rich rock (leach-

2318

M. 0. Siegel and S. Anderholm

103

2J

+

Ll

A

Lib3

10"

/ji

A =i 2

E f g a

10s -

-

10”

% 2

161_

?&

Lb

A *A

100

4 -+JT

10.’

? 0

+

A

+

+

*B

i;s m.i

A

+

+

10-s

-

f3 A

-f

10-6

y+

I

I

I

I

I

I

1

2

3

4

5

6

lo-’ 7

ionic Strength (mot&) FIG. 15.Relationship between Li concentrations (mg/L), Li/CI weight ratios ((mg/L)/(mg/L)). of Culebra groundwater samples. From SIEGELet al. ( 199 I )

limited) and the accessibility of silicates to groundwaters. Calculations by SIEGEL et al. ( 199 1) show that Culebra waters are probably all undersaturated with respect to amorphous silica. Saline waters in hydrochemical-facies Zone A lie close to quartz saturation; fresh waters from Zone 3 have silica contents between tridymite and amorphous silica saturation. No pattern is apparent for water samples from Zone C. The spatial pattern of silica concentration contours is different than patterns for elements which are assumed to be dissolved from halite and associated evaporite salts (cf. Fig. 4). In general, silica concentrations are highest in the west and decrease to the east. The lack of solubility control on silica concentration and this spatial trend suggest that the distribution of silica is strongly controlled by the rate of supply of silica to the ~oundwater (leach-limits ) . Source rocks in the western part of the site may contain more silica that is accessible to leaching than do rocks to the east. This may be due to the greater amount of evaporite dissolution and concomitant concentration of residual detrital silicates in the western and southern parts of the site compared to the eastern section. This hypothesis is supported by the studies of the mineralogy of samples from intact core described previously. Clays comprise about 3-5% by weight of the bulk mineralogy of intact Culebra core samples and are commonly associated with silt-sized quartz. Because silica will be released by dissolution of both quartz and clays, the distribution of silica as well as Li and B may be controlled by the accessibility of the clays to groundwaters. In areas where clays are abundant, both dissolution and ion exchange may ail&t water chemistry. The extent of release or uptake of Li and B by clays could depend on a large number of solution/substrate parameters including surface area, salinity, temperature, and past diagenetic history. Ion exchange data are not available for WIPP clays and brines; however, the info~ation provided above suggests that ion exchange probably occurs in Culebra waters and may exert some control on the distributions of Li and B in this aquifer.

Cornpieme~~~

and ionic strengths

Nature of Modeling Approaches

Each of the interpretative methods described in this paper has a unique set of advantages and limitations when evaluated in iight of potential errors in the available groundwater chemistry dataset (e.g., pH measurements) and uncertainties in the hydrological models for the WIPP site. Interpretative techniques such as Piper diagrams and principal component analysis ( PCA ) examine element ratios and correlations, and use the largest possible datasets. These techniques require no ad hoc assumptions about hydrologic llow fields or chemical equilibria and are relatively insensitive to errors in individual chemical analyses. Because of this lack of constraints, however, the results of the principal components analysis can only be used to suggest possible reactions; inde~ndent confirmation must be obtained for the proposed mechanisms. Other limitations of the PCA presented here are the small size of the dataset and the large relative uncertainties in the concentrations of some of the minor and trace solutes. In this work, calculations of mineral saturation indices and reaction path calculations were carried out with the PHRQPITZ code to evaluate the feasibility of geochemical processes suggested by the PCA. The relationships between K, Mg, Ca, SOa, and chloride in the Culebra waters were in quali~tive agreement with those generated by the PHRQPITZ simulations. The resulting reaction path model is constrained by thermodynamic relationships but not by individual chemical analyses or assumptions about the hydrologic flow field. Finally, models of net geochemical reactions along modern flow directions were calculated with the NETPATH code. These calculations required assumptions about the hydrologic flow field, were very sensitive to errors in individual chemical analyses, and are not constrained by thermodynamic or geologic data. All posaibIe chemical reactions consistent with the chemical analyses and proposed phases were calculated for several pairs of wells. Reaction sets that include dissolution

Geochemical evolution in a dolomite aquifer of halite, carbonate and evaporite salts, and ion exchange were consistent with the results of the PCA, the reaction path calculations with the PHRQPITZ code, and with observations of local mineralogy. The agreement among the results of the different types of calculations and the consistency of the proposed mechanisms with the available geological data lead to some assurance that the processes that have strongly influenced the chemical evolution of water in the Culebra dolomite have been identified.

SUMMARY

AND CONCLUSIONS

A model for the post-Pleistocene hydrochemical evolution of the Culebra Dolomite Member of the Rustler Formation has been formulated from several complementary geochemical modeling approaches. Salinities of water in the Culebra range roughly from 10,000 to 200,000 mg/L. The concentrations of the major solutes Cl, SO4, Mg, Ca, and K, as well as the Mg/Ca molar ratio, are correlated with the Na concentration. Three kinds of factors were extracted by R-mode principal component analysis from the major and minor solute data. The first factor (“salinity” factor) is dominated by Na, K, Mg, Br, and Cl. A second factor (“sulfate” factor) also has a strong Na-Cl influence, but is dominated by Ca, sulfate, bicarbonate alkalinity, and Sr. Geologic data from the site suggest that these two factors represent addition of solutes by the dissolution of halite, gypsumlanhydrite, and carbonates. A third factor (“silicate/bicarbonate” factor) showed the interelement correlations independent of effects attributable to halite dissolution. This factor contains two groups of correlated elements: Mg, bicarbonate alkalinity, and Si02 form one group, and pH, B, and Li form another group that is inversely correlated to the first group. This pattern of element associations might reflect sorption of B and Li by clays, release of Mg and SiOz by dissolution of authigenic clays, and carbonate diagenesis. Waters from the Culebra are saturated with respect to gypsum and undersaturated with respect to halite; calculated saturation indices for calcite and dolomite indicate saturation for some wells, but many of the waters show apparent supersaturation. It is likely that calcite and gypsum are dissolving in areas where halite dissolution occurs due to increases in solubilities with increasing ionic strength over the range of O-3 molal. The Mg/Ca ratios and calculated saturation indices of the waters are consistent with increasing degrees of dolomite supersaturation with ionic strength. The PHRQPITZ code was used to calculate the compositions of waters that would be produced by several hypothetical reaction paths driven by salt dissolution. The relationships between K, Mg, Ca, and SO4 with chloride are consistent with a partial-equilibrium model in which saturation is maintained with respect to calcite and gypsum while these solutes are added to the Culebra by dissolution of halite containing minor amounts of polyhalite and camallite. Recharge waters entering the Culebra may reach equilibrium with dolomite initially but as the ionic strength increases due to dissolution of halite and the accessory minerals, the solutions become increasingly supersaturated with respect to dolomite. The relationship between modem flow direction and groundwater compositions in the Culebra dolomite within

2319

the WIPP Site was examined using the NETPATH code. The mass transfers required to account for changes in the chemical compositions of groundwaters were inferred from analyses of waters sampled in a series of wells along postulated modem flow paths. In general, two kinds of models are consistent with the observed changes in water composition. One set of models requires ion exchange; the second set of models requires dissolution of relatively large amounts of evaporite minerals that have not been reported in the Culebra or are not common in the Rustler layers that stmtigraphically bound the Culebra. The abundance of clay minerals within and adjacent to the Culebra suggests that models involving ion exchange are reasonable. The work described in this paper illustrates the complementary use of several different modeling approaches. Each of these techniques has a unique set of advantages and limitations. Techniques such as Piper (trilinear) diagrams and principal component analysis (PCA) use the largest possible datasets, require no ad hoc assumptions about hydrologic flow fields or chemical equilibria, and are relatively insensitive to errors in individual chemical analyses. They can be used to suggest possible reactions which can then be evaluated using thermodynamic reaction path calculations and models of net geochemical reactions along modem flow paths. A unique model describing the chemical evolution of the groundwater system might not be obtained from any single technique; however, if it can be shown that the different geochemical models provide consistent results, and if the processes suggested by the model simulations are consistent with geohydrological data from the site, then we have some confidence that models are useful in identifying the mechanisms that have strongly influenced the chemical evolution of the groundwater. Acknowledgments-The authors are grateful to S. J. Lambert for allowing them to abstract his summaries of the geohydrology and isotope geochemistry of the WIPP site for this paper. K. Robinson’s careful reviews of the hydrochemical data produced a reliable data set for the calculations described in this paper. T. Sewards carried out the mineralogical studies that provide important constraints on the interpretation of the geochemical models. L. N. Plummer (USGS) allowed us to use an unreleased version of the PHRQPITZ code for the saturation index and reaction path calculations and provided some instruction in its use. Suggestions by B. F. Jones and L. N. Plummer (USGS) were very helpful in interpreting the results of the model calculations presented here. I. J. Hall (SNL) assisted with the principal component analysis. Comments from reviews of the draft manuscript by L. N. Plummer, B. F. Jones, R. Pabalon (SWRI), and J. H. Krumhansl (SNL) resulted in an improved final manuscript. This paper was prepared for the US Department of Energy under Contract DE-AC04-76DO0789. Editorial handling: H. Ohmoto

REFERENCES BACKW., HANSHAW B. B., PLUMMERL. N., RAHN P. H., RGHTMIRE C. T., and RUBIN M. ( 1983) Process and rate of dedolomitization: Mass transfer and 14Cdating in a regional carbonate aquifer. GSA Bull. 94, 1415-1429.

BARR G. E., MILLERW. B., and GONZALESD. D. ( 1983) Interim Report on the Modeling of the Regional Hydraulics of the Rustler Formation. Sandia National Laboratories Report SAND83-039 1.

BASSETT R. L. ( 1982) Geochemistry and mass transfer in brines from Permian Wolfcamp Carbonates in the Palo Duro, Dalhart,

M. D. Siegel and S. Anderholm

2320

and Anadarko Basins. In Groiogy

Groundwu/er

Euro Basin, Texas Panhandle:

Lab. Rep. SAND86-1054.

and Geohydrology o/‘the Palo A Report on the Progress of Nuclear

U’aste Isolation Feasibility Studies II 9SI j. Texas Bureau of Economic Geology Report. BATES R. G. ( 1975) pH scales for sea water. In The Nature of&awater; Dahlem Workshop Report (ed. E. D. GOLDBERG). pp. 3 13338. Verlagsgesellschaft. BATESR. G. and CULBERSONC. H. ( 1977) Hydrogen ions and the thermodynamic state of marine systems. In The Fate of‘ Fossil Fuel CO2 in the Oceans (ed. N. R. ANDERSON and A. MALAHOFF), pp. 45-6 I. Plenum. BEA~JHEIMR. L. ( 1987 ) Interpretations of Single- M’ell Hydraulic Tests Conducted at and near the Waste Isolation Pilot Plant (WIPP) Site, 1983-1987. Sandia Natl. Lab. Rep. 87-0039. BODINE M. W.. JR., JONES B. F.. and LAMBERT S. J. ( 1989) Nor-

mative

analysis

of groundwaters

of

from the rustler formation. In

Nr~drogc,ochcmicul Studies theRustler Formution Rocks in the WIPP .4rea, Southeastern New Me%)

und Relaled

(ed. M. D. SIEGEL. et al.), Chap. 4. pp. I-140. Sandia Natl. Lab. Rep. SAND880196. CHAPMANJ. B. ( 1988) Chemicaland Radiochemical C’haracteri.stic.s o/‘GronndH~atrr in the Cnlebra Dolomite, Southeastern .New Mesice. New Mexico Environmental Evaluation Group. Report EEG-

3’). COLLINS A.

G. ( 1975) Geochemistry of’Oil Field Water.s. Elsevier. DREV~;RJ. 1. ( 1982) The Geochemistry cfN&rrul A’utcn PrenticeHall. EVANSG. V., OTLET R. L.. DOWNINGR. A.. MONKHOLJSE R. A.. and RAEG. ( 1979) Some problems in the interpretation of isotope measurements in United Kingdom aquifers. Proc. ofthe Intemational Sympo.sium on Isotope Hydrolog!!; Intl. Afomir, Energ!. Awnc~~ STI/PUB 493. V2 (Vienna), oo. 679-706. GOI~BE~G S.‘and GLAUBIG‘R. A. ( i 9k6) Boron adsorption and silicon release by the clay minerals kaolinite. montmorillonite, and illite. Soil Sci. Sot. Amer. J. 50, 1442-1448. GONZALEZD. D. ( 1983) Groundwater Flow in the Rustler Formatron, M’u.ste I,wlution

Pilot

Plant

(U'IPP) ,

Southeast

Nevv Me.rico

(SSNM) Inrerim Report. Sandia Natl. Lab. Rep. SAND82-1012. HARDERH. ( 196 I ) Einbau van Bor in detritische Tonminerale. Experimente zur Erklarung des Borgehaltes toniger Sediments. Gcochim. C’osmochim. Acta 21, 284-294. HARRISSR. C. ( 1969) Boron regulation in the oceans. Natwe 223, 290-291. HARVILC. E. and WEAREJ. H. (1980) The prediction of mineral solubilities in natural waters: The Na-K-Mg-Ca-Cl-S04-Hz0 system from zero to high concentration at 25°C. Gcochim. c’osmochim. Actn 44, 98 l-997. HARVIEC. E., MOLLERN.. and WEAREJ. H. ( 1984) The prediction of mineral solubilities in natural waters: The Na-K-Mg-Ca-HCl-SO,-OH-HCOj-COj-COZ-HZO system to high ionic strengths at 20°C. Geochim. Cosmochim. Acta 48, 723-75 I. HITC‘HONB., BILLINGSG. K., and KLOVANJ. E. ( 197 I ) Geochemistry and origin of formation waters in the Western Canada Sedimentary Basin: III. Factors controlling chemical composition. Gcochim

Cosmochim.

Acta 35, 567-598.

HOLSERW. T. ( 1979) Trace Elements and Isotopes in Evaporites. In Marine Minerals (ed. R. G. BURNS), Vol. 6. Chap. 9, pp. 295346. Mineral. Sot. Amer. JONES B. F. and WEIR A. H. (1983) Clay minerals of Lake Abert. an alkaline, saline lake. Clap Clay Minerals 31, 16 l-1 72. KI OVEN J. E. ( 1975) R- and Q-Mode Factor Analysis. In Concepts in Gr~stutrstrc.s (ed. R. MCCAMMON ), pp. 2 l-69. Springer-Verlag. KRUMHANSLJ. H., KIMBALLK. M., and STEINC. L. ( 1991) A Review of M’IPP Reporitor~~ Clays and their Reiutionship to Clays of AdJawnt Struta. Sandia Natl. Lab. Rep. SAND90-0549.

LAMBERTS. J. ( 1978) Geochemistry of Delaware Basin groundwaters. In Geology and Mineral Deposits of‘ochoun Rocks in Delaware Basin und Adjuccnt .Ireas (ed. G. S. AUSTIN), pp. 33-38. N. M. Bur. Mines & Min. Res. Circ. 159. LAMBERTS. J. ( 1983) Dissolution of’Evaporlte.s in and Around the Delaware Basin, Southeastern New Me.xico and West Texas. Sandia Natl. Lab. Rep. SAND82-0461. LAMBERTS. J. ( 1987) Feasibility study: Applicability of‘Geochronolqprc. Methods

Involving

Radioearhon

und Othcsr Nrrclides to the

IIJ!dro/ogy

($/ho

Rwrler

E‘ormation.

Sandia

Natl.

LAMBERTS. J. ( I99 I ) Isotopic constraints on the Rustler and Dewey Lake groundwater systems. In Hydrogeochemical StudicJs of the Rustler Formation wstern New, Meriw

and Related

Rocks in the WIPP

Ircu. South-

(ed. M. D. SIEC;EI.et al.). pp. 5. I-S-79. Sandia Lab. Rep. SAND88-0 196. LAMBERTS. J. and CARTERJ. A. ( 1987) Cranium-l.co/ope Sy.stema/rc:r in Groundwa/er.s ofthe Ru.stler Formation, Basin. Sinrthca.~twn ,Ycw .Meyic,o. Sandia Natl.

Northern

Delawure

Lab. Rep. SAND87-

0388.

LAMBERTS. J. and HARV~VD. M. ( 1987)

StablcI.sotopr Gcochemi.strJ’ o~‘Ciroirndl~ut1,r.s in the D&ware Basin of Southcu.stcrn Neu Mwco. Sandia Natl. Lab. Reo. SAND87-0 138. LAMBERTS. J. and ROBINSONK..L. ( 1984) Field GeochemicalStudits of Groundn,ater.s rn Nu.sh DraLv, Southeustern .New .Ue.xico.

Sandia Natl. Lab. Rep. SAND83-1122. LAPPIN A. R. ( 1988) Summur~~ of Site-Characterizatwn Conducted From 1983 Through Plant (M ‘IPP) Site. Sorrtheustrrn

Studies Pilot

1987 at the Waste Isolution New Mexico.

Sandia Nat]. Lab.

Rep. SAND88-0157. LAPPIN A. R. and HUNTER R. L.. eds. (1989)

Systems ,4nul.vsi.s. Long-term Rudionuclide Trunsport, und Dose .4sse,ssment.s. Waste /solution Pi/o/ Plant I WIPPJ, .southeaa.stern Ne~v Me.rico. Sandia

Natl. Lab. Rep. SAND89-0462. LAVENIJ~ A. M.. HAUL A.. and KELLEY V. A. (1988) Numerical Simulation of Ground-water Flow in the Culebra Dolomite at the U’u.ste Isolution Pilot Plant ( WIPP) Site: Second Interim Report.

Sandia Nat]. Lab. Reo. SAND88-7002. LAVENUEA. M.. CA&MAN T. F.. and PICKENSJ. F. Gronndrwrcr C‘alibrution.

Flow Modeling

of the Culebra Dolomite:

(1990) Vol. 1: Model

Sandia Natl. Lab. Rep. SAND89-70681 I. LYONM. L. ( I989 ) Annual Water Quality Data Report for the Waste I~~lut~on Pi/o/ P/unt. Department of Energy Report DOE/ WIPP 89-00 I

MARIETTAM., BERTRAM-HOWERY S. G., ANDERSOND. R., BRINSPER K. F.. GUZOWSKIR. V.. IUZZOLINOH., and RECHARDR. P. ( I989 ) Perfimnunce .-l.~.se.ssment Methodology Demonstration. Methodology Development fi)r Evahlating Compliance With EPA 40 CFR 191. Subpart B.,/or the Wa.s/e Isolution Pilot Plant. Sandia

Natl. Lab. Rep. SAND89-2027. MEIJERA. J.. LOIXAMAJ. L., and PEARSONF. J. ( 1987 ) Modeling o/~Gronnd- Mhtw FIoM. rn the c’ulebra Dolomrte at the Waste Iso/u/ion Pilot Plunt ( M’IPP) Site: Interim Report; Appendix E (ed. A. HAAG et al.). Sandia Natl. Lab. Reo.. SAND86-7 167.

MERCERJ. W. ( 1483) Geoll!~drolog~~ Pilot

Plant

Site, 1,o.s Medanos

ofthi Proposed Waste Isolation A reu, Southeastern New Mexico.

USGS Water-Resources Investigation Report 83-4016. MERCER J. W. and ORR B. R. (1979) Interim Datu Report on the Geoh~~drolqq~ of IMP Proposed U’astc Isolation Pilot Plant Site. Southeast Nrtc, ,Me.wo USGS Water-Resources Investigations

Report 79-98. MERCERJ. W. and SNYDERR. P. (1990a) Drillholes

11-14 und I$- 13 / Waste Isolution

Sandia Nat]. Lab. Rep. SAND89-0202. MERCERJ. W. and SNYDERR. P. ( 1990b) Drillholes W’IPP).

at the H- 11 Comp1e.r

I Waste

Busic Dutu Report .fitr Pilot Plant- WIPP). Basic Datu Report .ltir Isolation Pilot Plunl-

Sandia Natl. Lab. Rep. SAND89-0200. MERCERJ. W. and SNYDERR. P. ( 1990~) Basic Drillhole

H-16 (Uhste

Isolation

Datu Report&r Pilot Plant- WIPP). Sandia Natl.

Lab. Rep. SAND89-0203. MERCERJ. W. and SNYDERR. P. (1990d) Basic Duta Report.for Drillholrs lr-I 7 and H- I8 ( U’uste Isolation Pilot Plant- WIPP). Sandia Natl. Lab. Rep. SAND89-0204. NOVAKC. F. ( 1992) An Evahrrrlion qf Radionuclide Batch Sorption Datu on C’ulebra Dolornitc$r .4queous Compositions ~hc Humun Inrrusion Scenurio,f~~r the Waste Isolution

Relevant to Pilot Plant.

Sandia Nat]. Lab. Reu. SAND86-09 17. PARKHURS~D. L.. THO~TENSON D. C., and PLUMMERL. N. ( 1980) PHREEQE-A

(bmputer

Program .fi,r Geochemicul

Calculutions.

USGS Water Resources Investigation Report 80-96. PEARSONF. J., JR., FISHER D. W., and PLUMMERL. N. (1978) Correction of groundwatcr chemistry and carbon isotope con-

Geochemical evolution in a dolomite aquifer sumption for effxts of CO, outsing.

Geochim. ~~srn~~~rn. AC&X

42, 1799-1807.

PEARSONF. J., JR., KELLEYV. A., and PICKENSJ. F. ( 1987) Preliminary Design for a Sorbing Trucer Test in the Culebra Dolomite at the H-3 Hydropad at the Waste Isolation Pilot Plant ( WIPP) Site. Sandia Nat]. Lab. Rep. SAND86-7 177, 4- 1I. P~TZERK. S. ( 1973) The~~ynamics of electrolytes. 1. Theoretical basis and general equations. J. Phys. Chem. 77,268-277. PITZER K. S. ( 1975) Thermodynamics of electrolytes. 5. Effects of higher-order electrostatic terms. J. Sol. Chem. 4, 249-265. PITZER K. S. and KIM J. J. ( 1974) Thermodynamics of electrolytes. 4. Activity and osmotic coefficients for mixed electrolytes. J. Amer. Chem. Sot. 96,5701-5707. PITZER K. S. and MAYORGAG. ( 1973) The~~ynamics of electrolytes. 2. Activity and osmotic coefficients for strong electrolytes with one or both ions univalent. J. Phys. Chem. 77,23OO-2308. PITZER K. S. and MAYORGAG. ( 1974) Thermodynamics of electrolytes. 3. Activity and osmotic coefficients of 2-2 electrolytes. J. Sol. Chem. 3,539-546. PLUMMERL. N. ( 1984) Geochemical modeling: A comparison of forward and inverse methods. Proc. First Canadian/American ConJ: Hydro1og.v(ed. B. HITCHONand E. I. WALLICK),pp. 14977. National Water Well Association. PLUMMERL. N. and SUNEQJISTE. T. ( 1982) Total individual ion activity coefficients of calcium and carbonate in seawater at 25°C and 35% salinity, and implications to the agreement between apparent and thermodynamic constants of calcite and aragonite. Geochim. Cosmochim. Acta 46,247-258. PI.UMMERL. N., PARKHURSTD. L., FLEMINGG. W., and DUNKLE S. A. ( 1988) A Computer Program incorporating Pitzer’s Equations for Ca~~ll~ation~fGe~hemical Reactions in Brines. USGS Water Resources Investigations Report 88-4 153. PLUMMERL. N., PRESTEMONE. C., and PARKHURSTD. L. ( 1991) An Interactive Code (NETPA TH),fir Modeling NET Geochemical Reactions Along a Flow PATH. USGS Water Resources Investigations Report 91-4078. POWERSD. W., LAMBERTS. J., SHAFFERS. E., HILL L. R., and WEART W. D., eds. ( 1978) Geological Characterization Report, Waste Isolation Pilot Plant ( WIPP) Site, Southeastern New Mexico. Sandia Natl. Lab. Rep. SAND78-1596, 2 vols. RAMEYD. S. ( 1985) Chemistry qf Rustler Fluids. New Mexico Environmental Evaluation Group Report EEG-3 I. RANDALLW. S., CRAWLEYM. E., and LYONM. L. ( 1988) Annual Water Quality Data Report for the Waste Isolation Piiol Plant. Department of Energy Report DOE-WIPP 88-006. REEDERR. J. ( 1981) Electron optical investigation of sedimentary dolomites. Contrib. Mineral. Petrol. 76, 148-157. REEDERR. J. and WENK H. R. ( 1979) Microstructures in low temperature dolomites. Geophys. Res. Let& 6, 77-80. ROBINSONK. L. ( 1992) Reference solute concentmtions in Culebm groundwaters at the Air intake Shaft (AIS) and Hydropads DOE1, H-3, H-l 1, H-14, H-15, H-17, and H-18. In An Evaluation of Radionuctide Batch Sorption Data on Culebra Dolomite for Aqueous Compositions Relevant to the Human Intrusion Scenario for the Waste Is(~~at~onPilot Plant (ed. C. F. NOVAK), Appendix B. Sandia Natl. Lab. Rep. SAND86-09 17. SEWARDST. ( 1991) Characterization ojFracture Surfaces in Dolomite Rock, Culebra Dolomite Member, Rustler Formation. Sandia Nat]. Lab. Rep. SAND90-7019. SEWARDST. ( 1993) Effect ofDisorder and Compositional Non-ldea!ity in Culebra Dolomite on the Gibbs Free Energy o~Formation. Contractor Letter Report to Sandia Nat]. Lab., Univ. New Mexico SAND90-2520.

2321

SEWARDST., WILLIAMSM. L., KEIL K., LAMBERTS. J., and STEIN C. L. ( 199 la) Mineralogy of the Culebra Dolomite. In Hydrogm chemical Studies of the Rustier Formation and Related Rocks in the WIPP Area, Southeastern New Mexico (ed. M. D. SIEGELet al.), Sandia Natl. Lab. Rep. SAND88-0196. SEWARDST., WILLIAMSM. L., and KEIL K. ( 1991b) Mineralogy of the Culebra Do~om~ieqf the Rusfler Formation. Sandia Natl. Lab. Rep. SAND89-7085. SEWARDST., GLENN R., and KEIL K. ( 1991~) Mineralogy of the Rustler Formation in the WIPP-19 Core. Sandia Natl. Lab. Rep. SAND87-7036. &WARDS T., BREARLEYA., GLENN R., MACKINNONI. D. R., and SIEGELM. D. ( 1992) Nature and Genesis of Clay Minerals ofthe Rustler Formation in the Vicinity qfthe Waste Isolation Pi/at Plant, Southeastern New Mexico. Sandia Nat]. Lab. Rep. SAND90-2569. SEYFRIEDW. E., JR., JANECKYD. R., and MOTTL M. L. (1984) Alteration of the oceanic crust: Implications for geochemical cycles of lithium and boron. Geochim. Cosmochim. Acta 48, 557-569. SIEGELM. D. ( 1992) Cal~i~~ions OfNet Ge~hemica~ M~s-Balance Reactions along Probable Hydrological Flow Paths in the Culebra Dolomite. Internal Sandia Natl. Lab. Lett. Rep. SIEGELM. D. and LAMBERTS. J. ( 1991) Summary of hydrogeochemical constraints on groundwater flow in the Rustler Formation. In H~~roge~hemical Studies ofthe Rustfer Formation and Related Rocks in the WlPP Area, Southeastern New Mexico fed. M. D. SIEGELet al.). pp. l-1-1-109. Sandia Nat]. Lab. Rep. SAND880196. SIEGELM. D., ROEINSONK. L., and MYERSJ. ( 1991) Solute relationships in groundwaters from the Culebra dolomite and related rocks in the WIPP area, southeastern New Mexico. In Hydrogeochemical Studies qfthe Rotter Formation and Related Rocks in the WlpP Area, Southeastern New Mexico fed. M. D. SIEGELet al.), pp. 2-1-2-164. Sandia Natl. Lab. Rep. SAND88-0196. SMITHC. L. ( 1976) Use of Lithium and Chloride Concentrations in Ground Water for Lithium Exploration. In Lithium Resources and Reqzi~rements bv the Year 2000 (ed, J. D. VINE): USGS Prqf: Paper f005, pp. 103-109. SMITH G. I. ( 1976) Origin of Lithium and Other Components in the Searles Lake Evaporites, California. In Lithium Resources and Requirements b.v the Year 2000 (ed. J. D. VINE); USGS Pro$ Paper 1005, pp. 92- 103. SNYDERR. P. ( 1985) ~isso~lltion of Halite and Gypsum, and Hydration ~fAnh~r~te io Gypsum, Rustler Formation, in the Vicinity ofthe Waste isolation Pilot Plant, Southeastern New Mexico. USGS Open File Report 85-229. SONNENFELDP. ( 1984) Brines and Evaporites. Academic. SPENCERR. J., EUGSTERH. P., and JONESB. F. ( 1985) Geochemistry of Great Salt Lake, Utah Ii: Plei~ocene-Hol~ene Evolution. Geechim. Cosmochim. Acta 49,739-747. Statistical Analysis System Institute ( 1982) SAS User’s Guide: Sfalislics. STEINC. L. and KRUMHANSLJ. L. ( 1986) Chemistry ofBrines in Salt from the Waste Isolation Pilot PIant ( WIPP), Southeastern New Me..uico:A Preliminary investigation. Sandia Natl. Lab. Rep. SAND85-0897. UHLANDD. W. and RANDALLW. S. ( 1986) Annual Water Quality Data Report. US Department of Energy Report DOE-WIPP-86006. UHLANDI). W., RANDALLW. S., and CARRASCOR. C. f 1987) Annual Water Qfia~ify Data Report, March, 1987. US Department of Energy Report DOE-WIPP-87-006.

7177 _. _-

M. D. Siegel and S. Anderholm

Table Al: Concentrationsof Solutest in Groundwatersfrom the Culebta Dolomite _~

Well(s)

Date2

Na (mUI-)

DOEI

@I/85 03/85 04/86 06/84 ous5 05i81 08/84 07l85 06/81 IO/81 08/‘85 05/a I 09l85 03/86 OH86 I Ii85 06185 08/85 02i86 03/86 02/87 08/80 02036 08i80 I Ii85 09l80 09l80 08/80 IV85 09/80 03/85

45800 18400 3570 17400 18000 6080 6150 5850 52400 52300 YIOO 18600 I8000 207 55.1 146 40400 49200 4360 28300 l9OcO 3160 3180 3620 4220 39200 15200 71400 94900 8570 200

DOE2 H-2A H-383 H-3B3 H4B H-K H4B H-5B H-SC H-SB HdB H&B

H-7Bl H-RB H-9B H-1183 H-12 P-14 P-17 WIPP-I3 w IPP-25 WIPP-25 WiPP-26 WIPP-26 W IPP-27 WIPP-28 WIPP-29 WIPF29 WIPP-30 fi8lC

SO.2. (mg/L) llO0 410 93.5 495 425 215 222 210 12% 1300 1350 450 375 7.0 3.83 6.85 943 1270 37.9 782 340 73.5 102 170 343 8060 485 IS600 23300 255 5.60

1610 I%0 167 829 783 455 505 427 2140 2150 2170 1080 1040 130 I57 137 1320 1980 840 1460

1730 I960 743 1550 1470 700 698 691 1710 1720 1700 2150 2040 587 548 5% 1700 1760 3520 1620

73600 34600 5310 29500 30300 7980 7950 7480 89500 89500 85400 33000 32300 320 305 194 65900 7!x!Rl 14500 48200 3lYxm 5250 6320 7200 8770 78500 24X00 138M)o 179000 14600 231

260 315 355 380 1900 555 5480 6500 460 152

905 1140 1240 1340 3210 II80 950 413 1140 588

Solute values have bee-n rounded as follows. - Na, K. Ca. Mg. Cl-SO, alkalinity values - all total dissolved solids (TDS) - allpHvalucs

7350 3950 2980 5130 4820 6230 5700 5520 7360 7560 7840 3980 3570 1850 1950 1900 7180 7210 I590 6020 4500 2500 2380 2480 2420 3830 43800 14aar 2OOco 4120 1990

Alkalinity3 as HCO; m )3

pH3 131200 6@lOrl 12900 55am 55800 21700 21200 20200 154400 154500 152600 59300 57400 3220 2830 3080 II7400 140500 24900 86500

7.1 7.0 8.0 7.4 7.4 8.0 7.8 7.7 7.9 7.9 7.4 7.0 6.9 7.3 7.3 7.4 7.2 7.2 6.8 7.5 6.6 6.9 7.2 6.9 7.1 6.4 6.5 6.1 5.9 8.8 7.4

45 67 57 52 71 75 69 80 86 50 96 94 120 96 110 54 53 I IO 64 -120 210 130 140 120 I20 210 160 40 110

l2lOa 13600 15100 17600 134700 46600 245400 324 IO0 29100 3270

56 34 5.6 29 26 42 48 43 62 64 49 34 34 0.57 0.085 0.24 47 76 72 72 39 2.6 3.4 3.2 3.9 28 73 45 61 IO 0.27

37 16 10 30 26 I8 20 I4 33 35 34 II IO 0.76 0.48 0.63 32 39 0.72 38 I1 I .5 1.7 1.4 1.6 2.3 5.8 4.4 5.2

0.64

0.47 0.22 0.53 0.40 0.39 0.49 0.40 0.77 0.77 0.81 0.44 0.45 0.10 0.12 0.18 0.62 1.2 0.28 0.87 0.35 0.20 0.22 0.24 0.23 0.33 0.30 0.78 0.70 0.27 0.17

26 38 9.5 23 30 I4 18 14 32 31 29 32 30 85 6.9 7.5 25 31 51 29 25 12 I7 I7 20 51 I6 29 I3 18 8.4

8.4 17 13 9.8 11 II 13 14 6.2 5.8 7.1 20 I8 47 29 27 7.2 30 8.5 34 33 33 35 23 36 22 IS 6.5 29

: 3 signilicant tigures : 2 sig. ligs. : lo nearest 100 mg/L or to 3 sig. figs. (for TDS
‘lk date on which samples were colkctcd. Alkalllity values and pH values were measured in lhc field when the samples were collected. Concentration of total dksolved sollds flDS). calculated by summlog the concentrations of majorsolutes. TableA3

TableA CompositionsofWatcrr Considc~dinM~s Tran.sferCalculations(inmgIL) Wells EkKiellt H-18 Ca Mg Na K Cl HC03 SO4

1,100 500

AIS 910 520

H-3 1,400 790

H-II 1,500 1,300

H-17 1,600 1,800

8,000 14,000 17,000 40,000 54,000 800

1,200

12,000 20,000 29,000 65.00

88,000

210

52 3,800

320

110 7,600

440

5, 4,800

54 6,400

50 7,200

Br

I4

27

29

50

76

B

16

30

26

32

43

7.63

74

7.25

70

PH Temp

755 19.6

16.6

22.5

23.3

21.3

I.02

1.04

I.04

1.08

I.10

Correlation MatrixforSoluteData'inC&bra Groundwaters

MkY

K

Na

Cl

SD4

B

Li

SiO2

Br

Sr

Ca

I 00

049

041

0 56

0.60

004

0 36

0 45

-025

0.63

0.94

-0.11

Mg

049

IO0

090

090

0.89

083

063

088

-0 58

0 85

071

-041

-0.15

K

041

090

100

094

093

0 84

0.66

0.80

-0.54

0.81

0 62

-035

-0.09

Na

0 56

0 90

094

100

099

081

079

0.88

-064

093

075

-022

-022

Cl

0 60

0 89

093

099

100

075

0.75

0 85

-0 60

094

0 78

-0 23

-0 19

SD4

004

083

0.84

081

0 75

100

0 75

0 83

-0 70

0 69

0 31

-0 I8

-0.22

B

"36

0 63

0 66

0 79

0 75

075

100

0 8b

-0.84

076

0 50

0.18

-0 57

0 85

083

0 86

100

-0 77

0.86

0 64

-0 I2

-039

-0 60

-0 70

-0 84

-077

I 00

-063

-0 39

-042

072

069

076

0 8b

-063

I 00

(1.81

-0 I5

-0 32 -0 26

L, SI02

045 -025

0 88 -0 58

080 -0.54

088 -0 64

(@ml)

PH

-022

Br

0 b3

0 85

081

093

0 94

Sr

0 94

071

062

075

078

031

0 50

064

-0 39

081

100

-020

-0 23

-0 18

0 I8

-0 I2

-044

-0 15

-0 20

100

-0.57

-0 I9

-022

-0 57

-0 39

072

-0 32

.O26

-0.57

too

PH HC03

-0 II -0 22

-041 -0 15

-0 35 -009

-022 -0 22

('C) Dens

HC03

Ca

l~he dataarethenatwallogarithms ofthesolute concentrations rnmg/L(except farpH)

Geochemicalevolutionin a dolomite aquifer

2323

Table A48 Utuotated &mode Factor loidndings Comm.

Element

Factor

IB

2B

3B

4B

5B

(hf)

1 Res. var’

Comm.’

Ca

2

1

Element

0.59

0.17

3

0.77

4

0.01

5

0.00

(h:)

0.97

Mg

092

0.25

-0.12

-0.14

-0.16

0.98

K

0.90

0.25

-0.2’

0.10

-0.14

0.95

Na

0.97

0.12

-0.03

0.13

0.01

0.99

Cl

0.96

0.16

0.01

0.15

0.00

0.98

SO4

0.81

El

0.84

Li SiO2 Br Sr PH HCO3

0.94 -0.74 0.94 0.77 -0.14 -0.37 61

-0.03 -0.38 -0.08 0.61

-0.56 -0.10 -0.12 0.11 0.12

0.059 0.19

0.58

-0.90 0.79 ‘8

0 ‘6

-0.03 0.04 -0.11 5.00 0.05 -0.03 0.34

-0.04 0.32 0.12 0.13 0.05 -0.04 -0.10

a.20

0.36

0.04

II

3

2

0.99 0.98 0.94 0.97 0.92 0.98 0.99 0.96

APPENDIX B. WATER COMPOSITIONS USED IN NETPATH CALCULATIONS Although groundwaters are electrically neutral, analytical errors and incomplete anaIyses frequently lead to an apparent charge imbalance in a water analysis (i.e., the sum of the charges of the cations do not equal the sum of the charges of the anions). The water compositions described in Table A2 were adjusted to eliminate this extraneous effect prior to calculations with NETPATH. This was accomplished by calculating the chemical speciation of the water using the computer program PHRQPITZ (Plummer et al., 1988; Version 0.2, Maclnnes Scale) and then calculating the total cationic and total anionic charges. The average charge (AC) was calculated as AC = total cation charge/2 + total anion charge/2; A charge balance factor was then calculated for anions and cations using the following formulae: Correction factor for anions = AC/total anion charge; Correction factor for cations = AC/totat cation charge. The concentrations of the charged species calculated with

Na

0.12

-0.10

0.76

-0.58

0.04

0.98

Cl

0.28

0.03

0.84

-0.36

-0.04

0.96

7%

Ci3

0.97

-0.01

0.14

-0.08

-0.0’

0.99

71%

MS

-0.03

0.23

-0.23

0.92

0.10

0.98

7%

K

-0.30

0.22

0.05

-0 ‘3

-0.25

0.98

9%

ST

0.94

-0.01

0.25

0.09

0.01

0.97

39%

PH

-0.01

-0.93

0.19

-0.26

0.07

0.98

91%

UC03

-0.1’

0.79

-0.04

0.03

-0.22

0.97

96%

-0.94

-0.18

.0.‘4

005

0.17

0.98

21%

Si02

0.18

0.93

0.04

0.09

-0.24

0.99

S8%

‘31

0.27

-0.16

0.86

0.08

0.13

0.96

‘9%

B

-0.19

-0.52

0.14

-0.55

0.47

0.97

43%

Li

-0.12

-0.34

0.07

0.07

0.90

0.99

18%

*,‘avar?

27

26

SO4

20

16

4%

‘1

‘Communality 2Variancc not correlated to tota’ dissolved solids (THIS); that is, the amount of variance remaining atk TDS was partialled out. 3Percent of Variance Explained by Each Factor

PHRQPITZ were multiplied by the correction factor to obtain the compositions used in the NETPATH calculation; therefore, the compositions of waters used in the NETPATH calculations differ slightly from those listed in Table A2. These corrections factors were small; for example, for the H-3 wet1 analysis, the correction factors for cations and anions were 1.0189 and 0.98 18, respectively. Species with concentrations of 10e5 mol/kg Hz0 or less were not included in the charge balance calculations. These included B species and several charged Ca, Mg, and S species. Concentrations of uncharged species of the major sofutes (Ca, Na, Cl, Mg, K, S) were less than IO-$ molfkg Hz0 and were not included in the NETPATH cakutations. Most of the total carbon in the water exists as the bicarbonate species; therefore, total carbon was corrected with the anion charge balance factor. The magnitude of the error introduced by this procedure is small and is a function of the concentrations of carbonic acid and carbonate ion. In the H-3 well, for example, the total carbon concen~tion is in error by 0:0007 mmol/kg H20 water because of the procedure used above.