WATER RESOURCES RESEARCH,

VOL. 36, NO. 5, PAGES 1277-1287, MAY 2000

Inverse chemical modeling and radiocarbon dating of palaeogroundwaters:The Tertiary Ledo-Paniselian aquifer in Flanders, Belgium W. J. M. van der Kemp Facultyof Earth Sciences,Free Universityof Amsterdam,Amsterdam Laboratoryfor Applied Geologyand Hydrogeology,Ghent University,Ghent, Belgium

C. A. J. Appelo2 Facultyof Earth Sciences,Free Universityof Amsterdam,Amsterdam

K. Walraevens Laboratoryfor Applied Geologyand Hydrogeology,Ghent University,Ghent, Belgium Fund for ScientificResearch,Flanders,Belgium

Abstract. Groundwatersamplesfrom the Ledo-Paniselianaquifer have been interpreted

for chemical reaction patterns, 14Cage,andrecharge conditions. Thisconfined Tertiary aquifer dipsNNE from its outcropin Belgium toward the North Sea over a length of-50

km.Conventional 14Cagesof thewatersamples rangefrom3 to over40 ka.Inverse chemical modeling wasdoneto correctthe•nCagesfor thechemical reactions in the aquifer,while accountingfor changesin the rechargewater qualityduringthe Holocene and late Pleistocene.The aquifershowsa zonal patternwith (goingupstream)Na-, K-, NH4-, Mg-, and Ca-HCO 3 water types.The pattern is a result of freshening:Ca displaces the salinecationsNa, K, NH4, and Mg from the aquifer'scationexchangecomplexin a

chromatographic sequence. Thelossof Ca2+ fromsolution bycationexchange isbyfar themostimportant reaction for dissolution of calcite, whichincreases theapparent 14C ageof thewatersamples. The14Cagefurthermore depends onopen/closed conditions of calcitedissolutionand CO2 gasexchangeand CO2 pressurein the rechargearea. It is

shown that•13CandCO2pressure in a soilareinterrelated andthatthechanges in CO2 pressurecan be includedin an inversemodel which considersvariationsin infiltration

waterquality. Theoverallcorrection for 14Cageis obtained byinverse modeling of water •3

qualityand •5 C, with optimizationon CO2 pressurein rechargewater usingPHREEQC [Parkhurst,1995]. The optimizedCO2 pressurefor the rechargearea varieswith age and is generallylower in the water sampleswith an age above13 ka. The lower CO2 pressureis

corroborated bylower;5180 valuesof thewater. 1.

an aquiferwhich infiltrated thousandsof yearsago.The reac-

Introduction

Inverse chemical modeling of groundwaterallows one to extractreactionpatternsfrom observedwater qualitiesin aquifers [Backet al., 1984;Plummeret al., 1990;Plummer,1992; Parkhurst,1997].The inversemodelsare basedon relating the concentrationsin a groundwatersample to an initial water qualityusinga user-specified set of reactions[Plummeret al., 1990, 1994; Parkhurst and Plummer, 1993; Parkhurst, 1997;

tionsin the aquiferare essential for correcting •4C agesof groundwatersamples. The modelspublishedso far [Plummeret al., 1983, 1990; Chapelleand Knobel,1985;Chapelleand McMahon, 1991;McMahon and Chapelle,1991;Aravenaet al., 1995] have not explicitlyconsideredhow the inversemodel parametersare affected by changes in recharge conditions. The climatic conditionsduring the late Pleistocenehave varied dramatically, giving rise to a potentiallylarge range of soil carbon dioxidepressuresand associated varyingcircumstances for calcite dissolution. It is obviousthat total inorganiccarbon(TIC) of groundwateris affectedby the soil conditions,but the fact

Aravenaet al., 1995; Chapelleand McMahon, 1991]. The extents of the variousreactionsin the set are obtainedby mole balancing [Parkhurstet al., 1982; Plummer et al., 1994; Parkhurst,1995]. The initial water quality may be that of a groundwatersampleupstreamalonga flow line in an aquifer, that{5•3C in soilwatercanalsochange withCO2pressure in or it maybe a guessat the water qualityin the rechargearea of the soil has not been evaluated in the inverse models. In fact,

analmost linearrelationexists between {5•3C of soilgasandthe •Nowat SubAtomicPhysics Department, Facultyof Physics and inverseof soil CO2 pressure[Cerlinget al., 1991]. This linear Astronomy,Utrecht University,Utrecht, Netherlands. relation allowsthe variationsin both soil CO2 pressureand 2Nowat Amsterdam,Netherlands.

{5•3C together to be included in theinverse model.

Copyright2000 by the AmericanGeophysicalUnion. Paper number 1999WR900357. 0043-1397/00/1999WR900357509.00

This paper considershow variationsin initial water quality can be included in inverse models. Each postulated initial water qualitywill lead to anothersetof reactionsin the aquifer. 1277

1278

VAN

-2.0

DER

KEMP

ET AL.'

CHEMICAL

MODELING

AND

RADIOCARBON

DATING

1/2

-3.0 1øg(Pcø2/pø) -3.3

I

,

I

,

I

,

aD= 13 D = M•3co2 + MairM:2co2 Mair_1 = 1.0044, 45 x29 44 +291 1/2

= 45+ 2944x 29

-10

where M stands for molecular

-15

(1)

mass.

The /513Cvalues of soil carbon dioxide are at least 4.4%0

enriched in •3C,buttheenrichment canbe more,depending

-20

-25

-30

0

500

1000

1500

2000

2500

3000

(1/ccO 2)' ppmV

on soil carbondioxidepressureand atmosphericcarbondioxide pressure,as was shownby Cerling[1984]with a diffusion model. The diffusiveflux of CO2 from the soil to the atmosphereover the concentrationgradientdc/dz can be written [D6rr and Mfinich, 1980]

Figure 1. The /513Cof soilcarbondioxideas a functionof

de j eKD • Iz=0

The set with the minimal residual error, calculatedby the inversemodel PHREEQC [Parkhurst,1995], is consideredto provide the optimal descriptionfor a givenwater sample.It

P(z) = Po exp (-Z/Zo),

= , (2) reciprocalsoil carbondioxideconcentrationin ppmV. Experimental data are from Dudziak and Halas [1996],Halas and Dudziak [1989],Cerlinget al. [1991],andRightmire[1978].The wheree is the air-filledporosity,K is a factorcorrectingfor the linesdrawnrepresentsoil-produced CO2with -29, -27, and tortuosity,andD is the diffusioncoefficientof CO2 in air. The -25%0. The commonpoint representsmodern atmospheric production of CO2 in the soilP(z) dependson depthz andis CO2. often representedby an exponentialequation[D6rr and Mt;innich, 1980; Hesterbergand Siegenthaler,1991; Dudziak and Halas, 1996]: (3)

whereP(z) decreases with increasingdepth(z is positivein a

willgiveanoptimal•4Cagecorrection andalsoyieldinforma- downwarddirection,and Po and Zo are constants).Following tion on the past rechargeconditionsbecauseCO2 pressure mustbe related to climaticvariations.An uncertaintystill resides in the choice of open versusclosed models for CO2 exchangeamongwater,calcite,andsoilgas.The idealcasesare illustratedand discussed in the paper. The conceptof inversemodelingwith variablerechargewater qualitieswasappliedto datafromthe Ledo-Paniselian (LP) aquifercollectedby Walraevens [1987].This confined,Tertiary aquifer extendsover a length of 50 km from its outcrop in Belgiumtowardthe North Sea.The water qualityin the aquifer showszonal bandsof Ca-, Mg-, NH4-, K-, and Na-HCO 3 waterswhich are a result of fresheningand cation exchange [Cardenaland Walraevens,1994; Walraevensand Cardenal,

Hesterberg and Siegenthaler [1991], the CO2 flux to the atmosphereand the CO2 production,at steadystate,canbe written

dj/dz + P(z) = 0.

With the exponentialP(z) functiongivenby (3), the differential equation(4) hasthe solution

C(Z) -- Cat m-{-C111-- exp (-Z/Zo)],

and late Pleistocene.

In this paperwe discussthe basictoolsthat were used,that

concentrationin ppmV. The isotopicmassbalancecan be written as

c(z)Rsoil--CatmRat mq- at)c•Rp[1- exp(-Z/Zo)],

1997].The reactiontransfers areusedto deduce:4C agesof the water samples.

(6)

whereR istheisotopic ratio:3C/•2Candthesubscript p refers to soil-produced CO2.This lastresultcanbe combinedwith (5) to obtain

is,therelationbetween/513C andCO2pressure in theinfiltrationareaandthereactions in theaquiferthataffect/5:3C and :4C ages.The reactions in the LP aquiferare reconstructed with the inversechemicalmodelPHREEQC [Parkhurst,1995,

(5)

where c• = (Po/e)(z•/,D) and½atm istheatmospheric CO2

1994].Conventional •4Cagesrangefrom3 to over40 ka and thusspanthe considerableclimaticvariationsof the Holocene

(4)

c(z) gsoil: Catmgatm -3-OID[C(Z ) -- Catm]gp.

(7)

Rearranging,introducingdelta notation, and writing Csoilc(z) leadsto

Catm(/513Catm q-1) q-aD(CsoilCatm)(/513% q-1)

/513Csoil --

c soil

2. Reactions Affecting•:3C and14C

(8)

in Groundwater

2.1. •13CandPcozoftheSoilAtmosphere

A test of (8), with experimentaldata from variousliterature sources,is shownin Figure 1. The lines drawn representCO2 13

It hasbeenshown byD6rrandMfinich[1980]that/513C of produced with /5 Cpvaluesof -29, -27, and -25%o. Note of the drawnlineswith soilCO2is enriched in 13Ccompared to the produced CO2 thatthe/5•3Cvaluesat the intercepts from degradationof organic matter or root respiration,be-

the/5•3C axisare4.4%0enriched compared to /5•3Cp. It is

cause•3CO2diffuses moreslowlyto theatmosphere thanthe evidentthat (8) offersan excellentapproximation for the carlighter12CO2. Theratioof thetwodiffusion coefficients (D) in bon 13 depletion as a function of reciprocal soil CO2 gas air leadsto a fractionationfactor (a) of 1.0044[Craig,1953],

concentration.

VAN

DER

KEMP

ET AL.:

CHEMICAL

MODELING

2.2. The •3C and Initial 14CActivityof Recharge TIC Threecasesof evolution of 8•3Cand14Cof recharge TIC are often distinguished in a sandyunsaturatedzone [Ingerson

RADIOCARBON

DATING

1279

Table 1. Atmospheric CO2Concentrations and813C Values During the Late Pleistoceneand Holocene As Derived

and Pearson, 1964; Pearson and Hanshaw, 1970; Deines et al.,

1974; Wigley, 1975; Mook, 1972, 1976; Fontes and Gamier, 1979;Reardonet al., 1980]. If the unsaturatedzone contains calcite,the chemicallyopensystemdissolutionof calciteis (1) full isotopicexchangebetweenTIC and soil CO2, denotedas

AND

From

Ice Cores

Period

Cco2, ppmV

Holocene Transition Glacial

279 243 198

8•3C,%0 -6.5 -6.6 -6.9

Averagevaluesare from Leuenberger et al. [1992] and Marino et al. (co,io),or (2) no isotopicexchange betweenTIC andsoilCO2, denotedas (co,ic). If the unsaturatedzone doesnot contain [1992]. calciteand the aquiferdoes,the chemically(and isotopically) closedsystemdissolutionof calcitein the aquiferis denotedas where AH2CO • istheactivity ofH2CO • inisotopic equilibrium (cc). with soil CO2. 2.1.1. Case (co,io). In this case,dissolvedtotal inorganic EmpiricalevidenceindicatesthatA rw< 100percentmodern carbon of infiltration water is determined by open system carbon(pmC); accordingto Vogel[1970],A rw - 85 +_5 pmC, equilibrationwith calciteand soil carbondioxide.The concenA rw = 65 +_5 pmC for recent trations of dissolvedcarbon dioxide, bicarbonate, and carbon- while Mook [1994] suggested prebombwaters.Furthermore,the useof one of theseempir-

atearefixed,whilethe8•3Cof thesespecies aredetermined by

icalA•wvaluesrequiresthatthe a13Cvaluesof the recharge

fractionationduring exchangewith the carbon dioxide from water mustbe corrected.This canbe doneby assumingpartial the soil atmosphere.After Deineset al. [1974] the equilibrium isotopicexchangebetween TIC and soil carbon dioxide as in distributionof carbon 13 among gaseousCO2 and aqueous the dissolutionexchangemodel by Mook [1972, 1976]. The

H2CO* 3, HCO•-, andCO32canbewrittenas Rco2(g)

initial14Cactivityof the recharge waterA rwcanbe approxiRco2(g)

matedby

Cl•CO2(g)-H2CO•(aq) -' RH2cO3•(aq) ' CI•CO2(g)-HCO•-(aq) = RI4co;(aq) ' Rco2(g)

(9)

A rw= (1 -- X•,'•rw _,x,• co,io co,ic + xArw ß

(14)

The choiceof A rw,at constantcarbondioxidepressure,fixes

thefactorx. Subsequently, •13Cof therecharged watercanbe

•CO2(g)-CO•-(aq) --Rco32_(aq),

approximatedfrom

where the a are fractionation factors, which have been re-

•j13Crw--(1 - x)•j13CCr•ø'iø q-x•j13CCr• ø'ic.

(15)

ported by Mook [1994] basedon data by Thodeet al. [1965], problem,withArwunknownand813Cof the Rubinsonand Clayton[1969],Emrich et al. [1970], Vogelet al. For the reverse [1970],andMook et al. [1974].Substitutionof the definitionof rechargewater known,the factorx can be calculatedfrom

a13Cin theformR = (a13C+ 1)RED Bleadsto thea13Cof the

a13Crw - a13CCr•ø'iø

aqueousspeciesi

•j13Ci(aq ) = aco2(g)_i(aq ) (•j13Cco2(g) q-1) -- 1,

X--'•13•.co,ic __•13•.co,io,

wherei represents aqueous H2CO•,HCO•-, andCO32-.813C of rechargeTIC is subsequently obtainedfrom

andArw canbe calculatedfrom (14). 2.2.3. Case(cc). In a chemically(and isotopically)closed

system, 813Cof TIC canbe writtenas

813CCOrw, iO= [C.2CO • 813 13 CracoW+ C•co•8 C•co•

+ cco•-a•3Cco•-]/c rw,

(16)

(10)

•j13('•cc / final initial,} •initiale13r-, -l•final •rw--' [ [Crw -- Crw /813Ccac03 q-Crw O !•,H2CO•J/Crw ,

(11)

(17) where Crw standsfor the TIC concentrationof the recharge where initial refersto equilibriumbetweenwater and soil carwater and the superscriptco,iostandsfor chemicallyopenand isotopically open.Theinitial14Cactivity of therecharge water bon dioxide(without calcite) and final refers to equilibrium betweendissolvedsoil carbondioxideand calcite(without a

A •ørw'iø canbeapproximated withanequation analogous to(11),

gaseous carbondioxidephase).The initial14Cactivityof the

usingfractionationconstantsa factorof 2 larger. rechargewater is simply 2.2.2. Case (co,ic). In the caseof a chemicallyopen but Ci•nitial isotopicallyclosedsystem,the concentrations of the chemical A • = Cfirv• speciesare identical to the previouscase. Since there is no nal,initial .cx TICß

(18)

isotopicexchange betweenTIC and soil CO2, 813Cof the HCO•-- and CO•--speciesin solutionis the averageof 2.3. RelationBetweenCO2Pressureand 8•3Cand •4C

a13CI_i2co • (asin(11))anda13Ccaco3; asa result, a13C ofTIC

of TIC

can be written

The aqueousconcentrationsof the carbon species,in the differentapproximations (co,io), (co,ic),and (cc), can be calculatedas a functionof CO2 pressurewith a chemicalspeciation computerprogram, for instance,PHREEQC, while the

as

813•.co,ic 13 - q-CCO32) vrw __[CH2CO•'8 Ca2co • q-•1(CHco• ø(•13CH2cO • q-•13Ccaco3)]/Crw.

(12)

a13Cof soilCO2canbeestimated with(8). Atmospheric CO2 concentrations and a13C values during the late Pleistocene and Assuming the14Cactivity for calciteto be zero,theinitial14C activityof the rechargewater becomes

rw• A.2CO•'[C.2co• + CCO•-)]/Crw, (13) Aco,ic q-•1(Cisco;

Holocene, required to do thesecalculations,are known from ice cores, [Leuenberger et al., 1992; Marino et al., 1992]. A summaryof thesedata is givenin Table 1. The possiblerange

1280

VAN

-2.0

DER

KEMP

ET AL.:

CHEMICAL

-3.0Iøg(PCO2/pø) -3.3 TIC(cc)

-3.5 ...........

......................• ...

• -10 -''•' -'";.••

MODELING

AND

RADIOCARBON

DATING

included in TIC, it follows from mass balance that TIC be-

comesgreatly enrichedupon methane production.Methane hasnot been reportedfor boreholesin the LP aquifer,and gas >' analysisfor a numberof water samplesrevealedno significant •'TIC(co,•o) amountsof methane.When complicationsdue to methanogenesiscan be neglected,dissolvedcarbonincreasesaccordingto

TIC(co,•c)•2(g)

(CH2O)106(NH3)16 + 212 H20

-15 -

-• 514H + + 106HCOj + 16 NH•- + 424e-.

(19)

-20 -

-25

Marine carbonatesshow a relative isotopicenrichmentin carbon13 of -1 to 2%o [Degens,1967]. Dissolutionwill add

...... 0

500

1000

1500

2000

2500

3000

3500

(1/cCo 2)' ppmV

this15•3C to groundwater TIC. Dissolution occurs whenacidis produced,for example,CO2 from oxidationof organicmatter,

for example, by cationexFigure2. The 15a3C of soilCO2andrecharge totalinorganic whenCa2+ is lostfromsolution,

carbon(TIC) as a functionof reciprocalsoil CO2 concentra- change, or when waters with different CO2 pressuresmix. tionfor soil-produced CO2with1513Cp = -27%0. Herecc Precipitationof carbonatesis possiblewhen acid is lost, for through reduction of SO42or Fe3+,orwhenprotons representschemicallyclosedsystemdissolutionof calcite;co,io example, representschemicallyand isotopicallyopen systemdissolution areexchanged, orwhenCa2+ or othersolidcarbonate forming of calcite; co,ic representschemicallyopen and isotopically ionssuchasFe2+ increase in concentration, for example, by closedsystemdissolutionof calcite. cationexchange or reduction in caseof Fe2+.Carbonisotope

fractionationbetweendissolvedbicarbonateandprecipitatesis approximately0.15%oat 10øC[Mook, 1994] and is thusinsigof 8•3Cof soilCO2 andcorresponding 8•3Cvaluesof recharge nificant.However,it canbecomesignificantif largeamountsof TIC, cases(co,io),(co,ic),and(cc),are shownin Figure2. TIC solid carbonatesrecrystallize,for instance,dissolutionof doin the fullyopensystemis 9.5%0heavier(at 11øC)thansoilgas lomite and subsequentprecipitationof calcite. The relative contribution of a reaction to calcite dissolution when CO2 pressureis smallerthan 0.01 atm becauseof fracor precipitationcan be estimatedin the first instancefrom rionation in the HCO•- species.At higher CO2 pressurethe

H2CO• species contributes increasingly to TIC, and15•3C of

eithertheprotons consumed or produced or theCa2+ gained

TIC becomeslighter. In the isotopicclosedsystem,TIC originateshalf from soilCO2 andhalf from calcite,with (for Figure

or lost in the reaction:

CaCO3 + H + Ca2++ HCO•.

2) a 15X3C of 1.67%o. The15•3C of TIC istherefore theaverage of the two contributions,while at higher CO2 pressurethe

There

will be additional

(20)

effects due to redistribution

of the

largercontribution of H2CO• lightens TIC again.Last,in the carbonatespeciesin responseto the pH changes.We givehere chemical andisotopic closed system the15•3C of TIC is domi- the most common reactions. nated by calcite at low CO2 pressure;while at higher initial 2.4.1. Loss of calcium. Calcium may be taken up in the CO2 pressure,it becomesgrossomodo the average of the exchangecomplex isotopicopen curveand of calcite,with againsomerounding 3 species. dueto increased importance of the H2CO* The corresponding initial•4Cactivities of recharge TIC in

the differentapproximations are shownin Figure3. The curves

Ca2++ 2X- -• CaX2

to be exchangedfor other base cations

for •4Ccanbe explained similarlyasfor •3C.TIC in the fully opensystem is2 pmCenriched in X4Cwithrespect to soilgas. In theisotopic closed system, x4Cistheaverage of soilgasand dead calcite, with an increasingcontribution of soil gas at higher CO2. In the closedsystem,dead carbonfrom calciteis the major contributorto TIC when the initial CO2 pressureis low, andA o is then small.

2.4. Reactions Involving•3C in the Aquifer The reactions whichcontribute to groundwater •j13Care oxidationor decompositionof organicmatter and reactionsof

(21)

(Na, K, 1/2Mg)X• X- + (Na+,K+, 1/2Mg2+). (22)

-3.0 1og(PCO2/pO) -3.3

110

I

-3.5

•IC(co,•o•

100

CO2(g) 90 80 70

• 60

TIC(co, •c)

4 50

solidcarbonates. Organicmatterhas15•3C = -25 +_5%oin a temperate climate. Microbial oxidation of organicmatter to

CO2via02, SO42-, or Fe(III) reduction willproduce CO2with approximately thesame15•3C astheorganic carbon. Methane producedfrom organiccarbon is stronglydepleted and has

15•3C ranging from-50 to -110%o.Methanefromthemarine environmentis usuallymore depletedthan methanefrom the freshwater environment,becausethe two possiblemethanogenic pathways,CO2 reduction and fermentationof organic matter, proceedat a different ratio. A reviewconsideringthis topic is given by Whiticaret al. [1986]. As methane is not

20

'''

•nC(cc)

10

0

0

500

!

!

1000

1500

2000

2500

3000

3500

(1/Cco2).ppmV

Figure3. Initial14Cactivity of soilCO2 andrecharge TIC as a function of reciprocal soil CO2 concentration for soil-

produced CO2with15•3Cp = -27%0. Corrections forfractionation were made based on an initial

14

C activity of the

atmosphere of 100%modernC (at 15•3C = -6.5%0).

VAN

DER

KEMP

ET AL.:

CHEMICAL

MODELING

The loss of Ca2+ from solution will induce dissolution of cal-

AND

RADIOCARBON

DATING

1281

ing to the discussed reactions(equations(19), (20), and (26)),

the•4Cactivity ofthesample withoutradioactive decay isgiven only10 to 20%of the Ca2+ loss,because pH increases when by calcitereleases CO32-to the solution. Observed calcitedissocite. Without proton exchangethe amountwhich dissolvesis

Crw

lution in aquifersin responseto Ca exchangecan be 100% of

And:• Arw.

the Ca2+ lossor of theNa+ gain[Chapelle andKnobel,1983; Hummeret al., 1994].The amountdependson the pH buffering capacityandcanbe calculatedfor a specificsolutionwith a forwardmodelsuchasPHREEQC [Parkhurst,1995].The dissolutionof calciteis largerthan calculatedfor a simple(Na, K, 1/2Mg)systembecausein the aquiferproton exchangeoccurs when the pH increases[Appelo,1994],

Ca2+q-2HX--> CaX2q-2H+.

(23)

For the caseof pureproton/Ca 2+ exchange the amountof calcitewhichdissolves canbe calculatedto be ashigh as 130 to 180% of Ca2+ loss.

2.4.2. Combined proton/redox reactions. Protons are consumedwhen iron-hydroxidedissolves,

FeOOH + 3H + -->Fe3+ + 2H20.

(24)

Protonsare producedwhen pyrite oxidizes,

FeS2+ 8H20--->Fe3++ 2SO]- + 16H+ + 11e-,

(25)

and when organicmatter oxidizesaccordingto reaction(19). The restriction

is that in the overall reaction

scheme the elec-

(29)

Notethatthe/•3C correction models, usuallyapplied,which have been reviewedby Fontesand Garnier [1979], formally cannotbe usedfor calculating And of the palaeowaterbecause these models assume a constant value for/•3C

of soil carbon

dioxide.

3.

Inverse Modeling

The chemicalreactionswhich determinewater quality and their relative importancecan be calculatedwith massbalance models or so-calledinversemodels [Parkhurstet al., 1982]. Originally,thesemodelswere limited in that masstransfersare constrainedby the stoichiometry of the reactions,becauseonly a set of linear equationsis solved.Parkhurst[1995, 1997] has recentlypublisheda new approachwhich allowsfor inaccuraciesin the analysisand for uncertaintyof the reactionscheme. The uncertaintyis addedasa separateterm (delta), andthe set of reactionequationsis solvedwith a Simplexschemeto the smallestvaluefor delta.This enablesone to (statistically)optimize the reactionschemeand to estimatethe optimalparam-

tronsare conserved, whichallowsoneto deducethe net proton eter values in the reaction scheme. A mass balance model which considers/•3C in addition to production.In the overallreaction,where sulfateS is reduced by organicmatter andprecipitatesin pyrite,sulfateis replaced variationsin CO2 pressurein the rechargeareawill enableone by bicarbonate: to deducethe amountof rechargedTIC, providedthat sufficient chemicaland isotopicdata are availableto use as con4 15 8 17 straints.We have used the statisticalinverse modeling ap106 FeOOH -+ • (CH20)•06(NH3)•6 + • SO•-+ 1-•H+ proachto optimizethe reactionschemein the Ledo-Paniselian 15 4 16 25 aquifer. ---> •-• HCO]+ • FeS2 q-•-• NH•-+ • H20. (26) The reactionshowsthat anionsare lost, togetherwith lossof

H +. Replacingsulfatewith the weakeracidbicarbonate will

4. Hydrogeologyand Hydrochemistry of the Ledo-Paniselian Aquifer

lead to further pH increase.Therefore this case of sulfate 4.1. Hydrogeology reductionby organicmatter can lead to calciteprecipitation. The givenreactions affect/5•3C of groundwater TIC, which The Ledo-Paniselianaquifer extendsfrom Flanders(Belfor an aquifer in which freshwaterand seawatermix, organic gium) to Sealand(Netherlands)over a distanceof-100 km. The aquifer is overlainby the Bartonian clay in the recharge matter oxidizes, and calcite dissolvescan be written as area and in the northernpart. The aquifer becomesphreatic •13G -- [ ( 1 -- fsw)C rw a13Crw q-fsw cswa13Csw south of the rechargearea. A map of the surveyarea, with + mcaco3•3Ccaco3 + mci_i2o•13Cci_i2o]/Cs (27) samplepointsfor hydrochemicaland isotopicanalyses[Walraevens,1987], is shownin Figure 4. A schematichydrogeowherefswis the seawaterfraction,rn the mass,transfer (mmol/ logical A-A' transect,indicated on Figure 4, from south to L), and Crwand csware the total inorganiccarbonconcentra- north through the recharge area southwestof Eeklo, is detionsof the rechargewater and seawater,respectively. Besides pictedin Figure5. At the northend of the sectionthe following the chemicalmassbalanceconstraints, (27) canbe usedas an stratigraphyis found:Quaternarysands,Boom (or Rupelian) additional 13C constraint. clay, Oligocene sands,Bartoon clay, Ledo-Paniseliansands, Paniselianclay,Ypresiansands,and finallythe Ypresianclay, 2.5. Radiocarbon Dating which is taken as the impermeablebottom layer of the conIf no TIC alteringprocesses takeplaceafter infiltration,then ceptualgroundwaterflow model. The Paniselianclay and the combination of the •4C activityafterinfiltration(no decay) Ypresian aquifer have an averagethicknessof 5 and 10 m, withtheobserved •4Cactivity (withradioactive decay)leadsto respectively.The thicknessesof the other aquifersand aquithe •4Cage,•4T, of a watersample: tardsvary alongthe transect. Two main rechargeareas exist;the most important one is And situated southwestof Eeklo at a topographicalhigh. Water infiltrates through the Bartoon clay, which locally varies in where the subscriptsnd an d standfor no decayand decay, thicknessbetween0 and 30 m, into the Ledo-Paniselianaquirespectively.However,when "dead carbon"dissolvesaccord- fer. The other recharge area is around Sint Niklaas where

14T = 8270 InAd'

(28)

1282

VAN DER KEMP ET AL.' CHEMICAL

MODELING

AND

RADIOCARBON

DATING

51ø13'N

51ø02'N

i

3o42'E

4o25'E

Figure4. Map of the surveyareawithsampling pointlocations for chemical andisotopicanalyses [Walrae yens,1987]and calculatedhydraulicheadsin the aquiferin its naturalprepumping state(before1920) [Walraevens, 1988].

observed in an upstreamdirectionin the Ledo-Paniselian aquifer [CardenalandWalraevens, 1994].The patterniscloselyakin to the sequenceobservedin the Aquia aquiferin the United 4.2. Cation Exchangeand Carbonate Reactions States[Chapelleand Knobel,1983;Appelo, 1994] and is the A sequencerangingfrom dilutedseawaterthroughsodium resultof cationexchangereactionsbetweenaquifersediments, (potassium) and magnesium to calciumbicarbonate wateris equilibratedwith seawater,and freshcalciumbicarbonatein-

rechargetakes place throughboth the Boom clay and the Bartoon clay.

S

N

m TAW 0-

Quaternary aquifer Boomclay

Oligoceneaquifer -100-

Bartoon clay

-200-

Ledo-Paniselian

aquifer

i

0

10 km

ß

Paniselianclay Ypresian aquifer Ypresian clay

Figure 5. Hydrogeological crosssectionA-A' of the aquifersystems andconfiningunits.

VAN

DER

KEMP

ET AL.'

CHEMICAL

MODELING

AND

RADIOCARBON

DATING

1283

SO42-

Table 2. Hydrochemicaland IsotopicAnalysesof Water Samples

Number

pH

TIC

Alk

22 44 84 20 83 11 2 26

8.07 7.90 8.38 8.02 8.23 8.59 8.34 7.48

11.7 10.2 11.0 10.6 12.6 12.9 10.2 8.30

11.5 10.0 11.1 10.5 12.6 13.2 10.3 7.69

28ph

8.31

6.20

Ca2+

6.19

Mg2+

Na+

K+

FeT

NH•-

C1-

0.22 0.44 0.19 0.31 0.12 0.07 0.25 1.61

0.48 0.94 0.28 0.42 0.18 0.11 0.32 1.17

23.9 44.7 13.2 13.0 24.9 15.0 38.9 1.96

0.70 1.05 0.56 0.69 0.45 0.38 0.63 0.69

0.0144 0.0081 0.0000 0.0039 0.0009 0.0018 0.0043 0.0104

0.0063 0.0422 0.0111 0.0242 0.0000 0.0000 0.0000 0.0632

18.1 34.3 5.37 3.91 12.4 3.23 28.2 0.62

0.20

•13C

14C

•180

0.81 1.90 0.31 0.18 0.34 0.49 0.44 0.11

-2.6 0.1 -5.2 -5.9 -2.5 - 1.8 -3.4 -13.6

2.5 0.5 4.5 1.8 0.8 3.0 1.2 45.2

-7.40 -7.03 -7.20 -7.00 -7.53 -7.06 -7.50 -5.57

0.00

- 12.0

42.0 16.6 33.6 39.9 12.4

-5.50

-6.56 -6.51 -6.71 -6.99

--8.1

37.7

--6.40

31 72 59 23

7.55 7.60 7.11 7.75

6.69 5.68 5.80 8.33

6.27 5.38 4.87 8.01

1.21 3.29 3.26 0.91

0.24

1.50 0.72 0.59 0.65

5.22 1.36 0.65 1.06 6.45

0.42

0.95 0.26 0.21 0.75

0.0029 0.0081 0.0183 0.0896 0.0335

0.0100

0.0742 0.0116 0.0269 0.0948

0.48 0.57 0.70 1.21 3.95

0.58 1.21 1.39 0.02

33ph

7.18

8.52

7.31

1.26

0.28

4.57

0.29

0.0450

0.1464

0.83

0.01

-

12.7 12.7 11.5 10.0

Source isWalraevens [1987].Chemical concentrations arein mmol/L.The•13Cand15•80 in percent milversus Peedee belemnite andstandard meanoceanwater,respectively. The •5x80 valuesarecorrected for seawater content(seawater, •5x80= 0%0).The X4Cactivities arein percent moderncarbon.Samplenumberswith subscriptph refer to phreaticwater samples.SubscriptT representstotal iron content.

filtrationwater.A chromatographic patterndevelops because pyrite oxidationin the infiltration area accordingto reaction the cations are displacedfrom the saline exchangerin the

orderof increasing selectivity fromNa+ firstto K + andNH•andfinallyMg2+ [Appelo, 1994].

(30) and subsequentreduction of sulfate in the aquifer by reaction(26). The reductionreactionwasthereforespecified

in theinversemodelfor samples witha SO]-/C1- ratiolower than seawater, while the oxidation reaction was indicated for

4.3. Inverse Chemical Modeling of Water Quality in the Ledo-Paniselian Aquifer

sampleswith a higherratio. The CO2 pressuresin the rechargearea were calculatedfor

Inverse chemicalmodelswere calculatedfor all samplesof Table 2 with PHREEQC-2 [Parkhurstand Appelo, 1999]. PHREEQC-2 includesisotopebalancesin additionto the uncertaintyestimatesin the massbalances.Mixing betweenfresh rechargewater (approximations (co,io), (co,ic),and (cc)) and seawaterwas calculatedfor reactions(19) to (26) and with

eachwatersample withPHREEQC-2usingstepsin logPco2

constraintsas indicated in Table 3. Furthermore, oxidation of

pyrite,which occursin the rechargearea in the coveringclay layer, was modeled by defining the "phase" (FeS2)4(O2)ls(H20)2.This givesriseto the reaction

(FeS2)4(O2)ls(H20)2---> 4Fe3++ 8SO]- + 4H+. (30)

of -0.1. Each CO2 pressuregivesrise to a numberof models, in our caseusuallyfrom 0 to 10 models.Each model is accompanied by a sum of the uncertaintiesgeneratedin the input concentrations (deltas[cf.Parkhurst,1995,1997]).It wasfound that the sumof the uncertaintiesfolloweda parabola,depending on soil CO2 pressure,and it wasjudgedthat the minimum in the curvewould correspondto the optimal model for the sample. An exampleof how the sum of uncertaintieschangeswith the logarithmof soil CO2 pressureof the infiltrationwater is shownin Figure 6 for sample31. There is a break-off edge associated with eachparabolawhichis indicatedas hatchedin Figure 6. Outside these edges, no models were found by PHREEQC-2 with the imposedlimits for the uncertainty.The

Subsequently,the ferric ions can precipitateas goethite.Reaction(30) was modeledas a completeoxidationreactionto separateit from the sulfatereductionreaction(26) whichmay take place downstreamin the aquifer. Furthermore, defining calculations were done for the three end-member cases of this phaseand reaction,insteadof a free 02 phase,prevents calcitedissolutionin the rechargearea for all water samples the generationof modelsin whichoxidationof organiccarbon from the aquifer. with oxygentakesplacein the aquifer.The direct oxidationof organicmatter in the root zoneis accountedfor by varyingthe 9 • CO2 pressureat infiltration. The presentlyavailable data do not permit separationof 8 CO, 7

Table 3. Mineral Phases,ExchangeSpecies,and

•6

Constraints

•5 O,

Chemical

Reactants

o 4 •3

Phases Exchangespecies

calcite,*goethite,pyrite,(FeS2)4(O2)•5(H20)2,* (cm20)106(Nm3) 16' CaX2, MgX2, NaX, KX, NH4X , and HX

Constraints

Na+, NH•-, K+, Ca2+,Mg2+,Fez,SO42-, S2-, CI-, TIC, Alk

2 1

0

-2

For 13C,seawater is 1.5_+0.1%o.For infiltration water,seeFigure

!

i

-1.9

-1.8

'

i

i

-1.7

-1.6

i

-1.5

i

i

-1.4

-1.3

i

-1.2

i

-1.1

-1

log(Pco 2/pO)

2 (error see Figure 2). For samplessee Table 2 (error _+ 0.1%o). Figure 6. Calculated sum of residualsas a function of log Calciteis 1.67 +_0.5%0;organicmatter is -25 _+2%0. *Phasesmay only enter the solution. pCO2 of the infiltrationwater for sample31.

1284

VAN

DER

KEMP

ET AL.: CHEMICAL

MODELING

AND

RADIOCARBON

DATING

aquifer. The ratio of infiltration volume to amountof calcite presentin the unsaturatedzone is such that equilibriumis

0.09

0.08

attained within

0.07

a few centimeters

of water flow in these soils.

Consideringthe rate of dissolutionof calciteunder thesecircumstances,open systemdissolutionis likely [Appeloand Postma,1993].Thereforeonly the modelresultsfor the chemical and isotopicopen systemare consideredin the following

0.06

0.05 0.04 0.03

discussions.

0.02

4.4.2. Contributions to groundwater TIC. The contributionsof rechargewater TIC, and of reactionsin the aquifer,to groundwaterTIC are shownin Figure 8. The contributionsare

0.01 0

-8

-7.5

-7

-6.5

-6

-5.5

-5

plottedagainst corrected •4Cageof thesample, whichis discussedin section4.5. The differencebetweenanalyzedground-

•80 / %0

waterTIC (indicatedbycirclesin Figure8) andthe sumof TIC Figure 7. The pCO2 at infiltration, derived from inverse in rechargewater and the dissolutionof calcitein the aquifer chemical modeling, asa functionof/•80 of thewater. (indicatedby crosses on a continuous line) is due to organic carbon oxidationin the aquifer (triggeredby sulfatereduction). Figure 8 showsthat TIC in the Pleistocenesamplesis 4.4. Results of the Inverse Model higher than in the Holocene waters. It also showsthat the 4.4.1. CO2 pressurein rechargewater. The optimal CO2 increaseof TIC is a resultof calcitedissolutionin the aquifer, pressures(thoseassociated with the lowestsumof uncertain- rather than in the rechargearea. The amount of calcitewhich ties)areplottedasa function of/•80 of thewater(corrected dissolvesin the rechargearea even declinesin the Pleistocene for seawater contentwith /•80 = 0%0)in Figure7 for the samplesbecauseof the lower CO2 pressurein the soil during three cases of calcite dissolution in the soil discussed. All the the colderclimaticconditionsin that period, aswasdiscussed models show that higher CO2 pressurescorrelate with less in section 1. The inversemodeldoesprovidethe relativecontributionsof negative/•80, implyingthat warmerconditions and higher organicproductivityare closelyassociated. It is of interestto the reactionswhich trigger the dissolutionof calcite in the thatis,lossof Ca2+ dueto cationexchange andproton note that the inversemodelingapproachhas thus provideda aquifer, methodto deducepalaeoclimatic conditionsby separatingthe productiondue to redoxreactions.The two typesof reactions reactionswhich have contributedto the groundwatercompo- are plotted cumulativelyversusamountof dissolvedcalcitein sition.Figure 7 showsthat a chemicallyclosedsystemfor the the aquifer in Figure 9. Figure 9 illustratesthat calcitedissocomplex rechargewater (casecc) has a muchhigherinitial CO2 pres- lutionis dueto theuptakeof Ca2+ in the exchange sure than the two other cases. The difference between the (reaction (21)), while the combinedproton/redoxreactions isotopicallyopenand closedcasesis minor,whichis due to the (reactions(19), (24), and(25)) addlittle or nothing.The small smalldifference in/•3C of TIC in recharge waterforthesetwo contributionof the redox reactionsin the LP aquifer also casesat the model-derived CO2pressures (compareFigure2). followsfrom the smallcontributionof the organiccarbonoxThe chemicallyclosedsystemrequiresa CO2 pressurein the idation reactionin the majorityof the groundwatersamples. It is rather conspicuous in Figure 9 that the net amountof rechargeareawhichis 5 to 10 timeshigherthan for the open case and appearsto be too high for the conditionsat this Ca2+lossishigher thantheamount of calcite whichdissolves.

o Sample TIC

I (a)+AquiferCC /

fla

•--

[]

•

,,,

[]

,,

13' ill'"

1 1

1

0

i

i

i

i

i

i

5

10

15

20

25

30

35

Corrected 14C age/ ka BP

Figure 8. Comparisonof contributionsof rechargeTIC and aquifercalcitedissolutionto sampleTIC as a functionof the correctedradiocarbonage. The differencebetweensampleTIC and the cumulativecontribution is due to organiccarbonoxidation.

VAN

DER

KEMP

ET AL.: CHEMICAL

MODELING

AND

RADIOCARBON

DATING

1285

14

- trl- .(a) + redoxreactions

o

10

= ...9.

8

'=

6

, -CaX 2 (a)

.g 4 E

2

i

0

I

i

i

i

i

i

i

i

i

1

2

3

4

5

6

7

8

9

10

CaCO3 dissolved / mmol.l 4

Figure 9. Cumulativeplot of the contributionsof cationexchangeand of redox reactionsto calcitedissolution in the Ledo-Paniselianaquifer. The amount of the cation exchangereaction is indicatedby the

diamonds(CaX2). The redoxreactions,(19), (24), and (25) in the text, are addedtogetherto the cation exchangereactionto yield the total reactionindicatedby the squares.The solidline showsthe l:t relation.

Theeffectisrelatedto thepH increase whichaccompanies the water agesare foundfor the caseof full chemicaland isotopic dissolutionof calciteand whichrequiresthat more calciumis

exchangeof rechargeTIC with soil CO2, while minimumages are obtainedfor the closedsystem.Ages derivedfrom chemamount removed isconsistently aboutt mmol/Lhigherforthe ically closedsystemdissolutionare comparableto those obfull rangeof calcitedissolutionfrom 2 to t0 mmol/L. The pH tained in the chemicalopen casewithout isotopicexchange. increasetriggersreleaseof protonsfrom the exchangecom- The correctionon apparentagecanbe ashigh as 15,000years, plex, and this accountsgenerallyfor almosthalf of the counter as is illustratedin Figure t0 for the fully open case. The exchange againstCa2+, the otherhalfbeingcaptured by the chemicallyclosedsystemis lesslikely for this aquifer as disremoved

from

solution

than can dissolve from

calcite.

The

base cations.

cussed before.

All latePleistocene watersin thisstudyhavemeasured •4C

4.5. Radiocarbon Ages

Conventionalradiocarbondating (no corrections)of the samplesshownin Table 2 leads to a maximum age >39 ka [Walraevens,t990a, b]. The optimized inversemodels,with

activitieslessthan 5 pmC (Figure t0), and thesesampleshave

the largestcorrections on •4C age.Hydrology, chloridecontent, and the extent of magnesiumexchangecan give addi-

recharge TIC (Crw)assource for •4Candaquiferreactions as tional cluesfor agedeterminationof samplesoriginatingfrom sourcefor dead carbononly, were used to obtain corrected the samerechargearea. The samplesare not generallyon the radiocarbonageswith (28) and (29). The maximumground- same flow line, but the correctedagesfor the open system 50

50

Apparent •4Cage 45

Corrected 14C age 40

•

oA

40

ß

35

m 30 •

30

25