Attribution: You must attribute the work in the manner specified by the author or licensor (but no in any way that suggests that they endorse you or your use of the work).Noncommercial ‒ you may not use this work for commercial purpose.No Derivative works ‒ You may not alter, transform, or build upon this work.Sharing ‒ Individual scientists are hereby granted permission, without fees or further requests to GSA, to use a single figure, a single table, and/or a brief paragraph of text in other subsequent works and to make unlimited photocopies of items in this journal for noncommercial use in classrooms to further education and science.


Modeling reactive geochemical transport in the subsurface is a powerful tool for understanding and interpreting geochemical processes in aquifer systems. Different conceptual models can include different combinations of geochemical processes. A limitation of current inverse models is that they are based only on one conceptual model, which may lead to statistical bias and underestimation of uncertainty. We present a stepwise inverse modeling methodology that can include any number of conceptual models and thus consider alternate combinations of processes, and it can provide a quantitative basis for selecting the best among them. We applied the inverse methodology to modeling the geochemical evolution in the Aquia aquifer (Maryland, USA) over 105 yr. The inverse model accounts for aqueous complexation, acid-base and redox reactions, cation exchange, proton surface complexation, and mineral dissolution and precipitation; identifies relevant geochemical processes; and estimates key reactive transport parameters from available hydrogeochemical data. Inverse modeling provides optimum estimates of transmissivities, leakage rates, dispersivities, cation exchange capacity (CEC), cation selectivities, and initial and boundary concentrations of selected chemical components. Inverse modeling with multiple conceptual models helps to identify the most likely physical and chemical processes in the paleohydrology and paleogeochemistry of the Aquia aquifer. Identification criteria derived from information theory are used to select the best among ten candidate conceptual models. In the final model, both proton surface complexation and methane oxidation are identified as relevant geochemical processes.


This study focuses on modeling the geochemical processes during the displacement of seawater by freshwater in the Aquia aquifer, a coastal-plain aquifer located in southern Mary-land. In his pioneer work, Back (1966) proposed the concept of hydrochemical facies to analyze the geochemical evolution and groundwater flow patterns in Atlantic coastal-plain aquifers. Freshening in coastal aquifer systems is accompanied by sequential elution of seawater cations from the sediment's exchange complex (Stuyfzand, 1993; Appelo, 1994). Flow conditions, transport parameters, and geochemical processes, such as aqueous complexation, acid-base and redox reactions, cation exchange, proton surface complexation and mineral dissolution and precipitation, have significant influence on the paleohydrology and paleogeochemistry of the Aquia aquifer. The complexity of the geochemical reaction patterns makes these processes difficult to model. Early attempts to interpret the chemical evolution of these systems used the mass-balance approach (Plummer and Back, 1980; Chapelle and Drummond, 1983), which was later implemented in computer codes, such as NETPATH, for modeling net geochemical reactions along a flow path (Plummer et al., 1994). Reactive transport codes, such as PHREEQC, allow modeling groundwater flow, solute transport, and chemical reactions simultaneously. Appelo (1994) reported a one-dimensional model for the Aquia aquifer that was calibrated by the trial-and-error method. In this paper, we show that inverse modeling of water flow and chemical reactive transport provides a methodology that allows us to select the best conceptual model for the interpretation of the geochemical evolution of the Aquia aquifer.

In the past decade, coupled models have been developed accounting for complex hydrological and chemical processes, and these have been applied to many case studies (e.g., Plummer et al., 1988, 1994; Yeh and Triparthi, 1989; Yeh, 2000; Engesgaard and Kipp, 1992; Appelo and Postma, 1993; Samper and Ayora, 1994; Parkhurst, 1995; Kohler et al., 1996; Keating and Bahr, 1998; Xu et al., 1999; Robinson et al., 2000; Sun et al., 2001; Zhu and Anderson, 2002; Gandhi et al., 2002; Saaltink et al., 2003; Lichtner et al., 2004; Molinero et al., 2004). In most cases, the chemical processes and parameters included in these models were calibrated by the trial-and-error method.

A mathematical and numerical methodology for solving the inverse problem of water flow, heat transport, and multicomponent reactive solute transport in variably saturated media was developed by Dai and Samper (1999a), verified by Dai (2000), and used successfully by Dai and Samper (1999b, 2004), Samper et al. (1999), Huertas et al. (2000), and García-Gutiérrez et al. (2001). A limitation of many current inverse models is that they are based on only one conceptual model, which may lead to statistical bias and underestimation of uncertainty (Neuman, 2003; Dai et al., 2004; Ye et al., 2004). We present here a stepwise inversion method, developed from Dai and Samper (2004), that can account for alternative conceptual-geochemical models of subsurface reactive transport and can provide a quantitative basis for selecting among them, while estimating flow and reactive transport parameters. The stepwise inversion strategy involves designing a list of the scenarios (or steps) that includes the most plausible combinations of geochemical processes, solving scenarios in a stepwise manner, and selecting the scenario that provides the best solution to the inverse problem. For example, the way to identify whether calcite dissolution/precipitation is relevant or not consists of solving the inverse problem under two alternate scenarios: (1) considering a geochemical system in which calcite is present, and (2) considering a geochemical system without calcite. After solving the two scenarios, it is usually possible to select the better result as the solution of the inverse problems and conclude whether calcite dissolution/precipitation is relevant or not. This stepwise strategy allows us to adopt a parsimony principle by testing models of increasing complexity and thus avoid overparameterization. It also helps us to identify the relevance of a given chemical process by solving the inverse problem under alternative scenarios in which the process is either occurring or not. There is no guarantee of a unique solution, but we present an example in which the stepwise inversion does lead to a parsimonious and tenable model for a complex geochemical system—a model that fits observations better than prior models.

The following section presents a model for the geochemical evolution of the Aquia aquifer system. The stepwise inversion strategy is then applied to the selection of the best conceptual geochemical model.


The Aquia aquifer is a coastal-plain aquifer located in southern Maryland (Fig. 1). It consists of a medium- to coarse-grained, medium- to well-sorted quartz sand (50%–75% volume fraction), glauconite, a hydrous iron-bearing aluminum silicate (20%–40%), and carbonate shell debris (1%–5%) (Chapelle and Drummond, 1983).

The Aquia Formation has been interpreted as having been deposited in a regressive marine environment. Figure 1 shows an outline of the Aquia aquifer with an estimated prepumping head distribution, which suggests a confined aquifer in the upstream part and a gradual loss of water in the downstream part of the aquifer (Chapelle and Drummond, 1983). Figure 2 shows a schematic cross section along a flow path shown in Figure 1. The transmissivity of Aquia ranges from 45 m2/d at outcrop areas to 90 m2/d at downgradient areas.

The Aquia Formation unconformably overlies the Upper Cretaceous Severn Formation, which according to Chapelle and Drummond (1983) can be considered as a confining bed. The Aquia Formation is overlain by the Marlboro Clay and the lower Nanjemoy Formation. The Marlboro Clay is a plastic clay composed predominantly of kaolinite and mixed-layer clays, while the Nanjemoy Formation is a glauconitic silty clay. The Piney Point Formation conformably overlies the Nanjemoy Formation and is a medium to coarse glauconitic quartz sand. This formation was truncated in upper areas by erosion that occurred during the withdrawal of the sea (Fig. 2). Upward leakage from the Aquia aquifer to the Piney Point Formation probably occurs via Pleistocene channels that cut through the confining beds (Chapelle and Drummond, 1983).

Hydrochemical data were collected from available water wells (Chapelle and Knobel, 1983). Near the outcrop of the Aquia aquifer, the water is of a calcium and magnesium bicarbonate type with low (6.5–7.5) pH. Downgradient from the outcrop, the water changes to a sodium bicarbonate type with higher (7.5–8.5) pH. The sodium-bicarbonate waters in the coastal-plain aquifers of the eastern United States are related to freshening of aquifers (Chapelle and Knobel, 1983).

Reported analyses of dissolved sulfate and CH4 indicate that their concentrations are very low (Chapelle and Knobel, 1985). Fe(II) concentrations are not reported, but most likely are too small to analyze (Plummer et al., 1994). Groundwater initially contains dissolved oxygen but becomes anoxic downgradient in the presence of lignitic material (Chapelle and Knobel, 1985).

Chapelle and Knobel (1983) presented the patterns of major cations along flow lines in the Aquia aquifer, and Chapelle and Drummond (1983) interpreted the observed geochemical patterns by using a mole-balancing model based on cation exchange and calcite dissolution and precipitation. They reported a cation exchange capacity (CEC) of 20 meq/100 g for glauconitic samples, which amounts to a CEC of 4–8 meq/100 g of aquifer sample. Shell material in the Aquia consists of calcium carbonate and magnesium calcite, Ca(1– x)MgxCO3(s), with x usually ranging from 0.01 to 0.18. Calcite cementation is abundant near the outcrop area and decreases downgradient.

Water quality in the Aquia aquifer was simulated by Plummer et al. (1994) using NETPATH and by Appelo (1994) using a one-dimensional model with PHREEQM. Appelo (1994) adopted a conceptual geochemical model based on that of Chapelle and Knobel (1983) and assumed that cation exchange and dissolution/precipitation of calcite are the major geochemical processes. He estimated solute transport and exchange parameters as well as initial and boundary concentrations with a trial-and-error method. He also proposed two possible mechanisms for the low pH measured in upgradient digestion and areas of the Aquia: methane/CO2 proton exchange. According to Appelo (1994) and Chapelle and Knobel (1985), the currently observed chromatographic pattern in the aquifer has been established during the last 105 yr. This is consistent with 14C activities in groundwater of the confined part of the Aquia aquifer, which suggest that most of this water infiltrated at least 30,000 yr ago (Aeschbach-Hertig et al., 2002).

Water quality in the aquifer shows zonal bands that approximately mimic hydraulic head contour lines. Figure 1 shows the average groundwater flow path. Available chemical data are projected onto this flow path for the purpose of their comparison to model results (see following).


The origin of sodium bicarbonate groundwater is generally attributed to the dissolution of calcium carbonate and cation exchange of Ca for Na on marine clays. It has long been recognized that a source of CO2 is required to react with calcium carbonate to form bicarbonate. Otherwise dissolution and exchange reactions would produce sodium-carbonate water with pH values above 10.3. Aquia sodium bicarbonate groundwater typically has pH values around 8.5. The source of CO2 has generally been attributed to organic material oxidized by microbial activity and has been a subject of considerable research (Chapelle and Knobel, 1983; Plummer et al., 1994). There is the possibility that sodium bicarbonate waters were formed by inorganic reactions involving exchange of sodium and protons for calcium.

Chapelle and Knobel (1983) noted that the sodium bicarbonate waters only form in the deeper, confined anoxic part of the Aquia aquifer. Therefore, aerobic oxidation of organic matter is unlikely, because there is no plausible source of O2 between the initial and final wells. So other electron acceptors must be considered. Plummer et al. (1994) considered and tested, using NETPATH, several conceptual models for the CO2 source, including: (1) aerobic oxidation of organic matter, (2) oxidation of organic matter by Fe(III) (goethite) reduction, (3) sulfate reduction, and (4) methanogenesis. Since Fe(II) concentrations are very low in ground-water, precipitation of magnetite or cation exchange of Fe(II) for Na must be postulated. In the absence of dissolved iron and H2S, it is necessary to consider sources and sinks of sulfur and iron, such as pyrite for iron and gypsum for sulfur. Although gypsum is not present in the sediments, pore waters in the adjacent confining units have been shown to be calcium-sulfate waters near gypsum saturation. Reduction of ferric hydroxide with Fe(II)-replacing Na on the exchanger requires an unrealistically large source of CO2 entering the groundwater (Plummer et al., 1994). Since there is no evidence of a reduced iron phase in the Aquia Formation, iron reduction alone is unlikely. Production of iron sulfide minerals during sulfate reduction could be a plausible mechanism. For sulfate reduction to occur in the Aquia aquifer, there must be a source and sink for sulfur and iron. FeOOH could be the source of iron, while the source of sulfate remains unknown because there are no evaporate beds or nodules of gypsum in the formation. The low CH4 concentrations in the Aquia aquifer eliminate the possibility of meth-anogenesis as a viable reaction. Chapelle and McMahon (1991a, 1991b) also concluded that observed methane production was inadequate to account for the observed CO2 production in the Black Creek coastal-plain aquifer. Plummer et al. (1994) concluded that Aquia groundwater may actually be explained rather simply without redox reactions. Proton exchange yields additional calcite dissolution, and the protons released from the exchanger eliminate the need for an additional CO2 source in the aquifer.

According to Chapelle and McMahon (1991a), a primary source of dissolved inorganic carbon in the Black Creek aquifer of South Carolina is CO2 produced by microbially mediated oxidation of sedimentary organic matter. However, groundwater chemistry data indicate that the available mass of inorganic electron acceptors is insufficient. Although sulfate concentrations are low in aquifer water throughout the flow system, sulfate concentrations are greater in the confining-bed pore water. From calculations based on Fick's Law, they concluded that possible rates of sulfate diffusion to aquifer sediments are sufficient to explain observed rates of CO2 production, thus eliminating the apparent electron-acceptor deficit.

Chapelle and Knobel (1985) assumed the dissolution of Mg-calcite. It is known that calcites of only low Mg content and aragonite are present in the clays, silts, and sands of the Paleocene-age Aquia Formation. Plummer et al. (1994) recognized the need to account for mineral sources of Mg such as low Mg-calcite and Mg-montmorillonite, and Mg-Na exchange. However, they ignored reactions involving Mg.


Our model represents a flow tube 60 miles long with a cross-sectional area of 1000 m2. Recharge takes place at the upgradient boundary, which coincides with x = 0. Similar to Appelo (1994), it is assumed that leakage into the confining layers is evenly distributed over the downstream half of the flow tube. This is a sensible assumption, because the first half of the aquifer is either a water-table aquifer or has weak confining conditions (see Fig. 2).

We use a two-dimensional code, INVERSE-CORE (Dai and Samper, 2004), to simulate flow through the flow tube per unit width. The flow domain is discretized by means of 70 triangular elements and 72 nodes. The nodal space along the streamline is 1.4 km. With the dispersivity in the aquifer set at ∼3.2 km (Appelo, 1994), the grid Peclet number is 1.1 and therefore acceptably small (smaller than 2). Similarly, the time increment is 3.5 yr, and with groundwater velocities of ∼0.6 m/yr, the Courant number is acceptable, since it is much smaller than 1. Two parameter zones are defined for transmissivity and dispersivity. Zone 1 includes the upstream part of the aquifer, while zone 2 corresponds to the leaky aquifer. Porosity is fixed and equal to 0.3 in both parameter zones. Water flow is assumed to be at steady state. Molecular diffusion within the aquifer is assumed to be several orders of magnitude smaller than hydrodynamic dispersion and therefore is negligible.

The Gaines-Thomas convention is used for cation exchange (Appelo and Postma, 1993). Na+ is used as the reference cation, and, therefore, its selectivity is equal to 1. Selectivities of other cations, such as K, Ca, and Mg (KNa/K, KNa/Ca, and KNa/Mg), are defined with respect to sodium and are estimated during the process of solving the inverse problem. Other chemical reactions include aqueous complexation and calcite and Mg-calcite dissolution and precipitation. The following eight species are taken as primary aqueous species: H2O, Na+, K+, Ca2+, Mg2+, Cl, HCO3, and H+. Relevant aqueous complexes were identified from speciation runs performed with EQ3NR (Wolery, 1983). They include: OH, CO32−, CO2(aq), CaHCO3+, MgHCO3+, and CaCO3(aq), MgCO3(aq), NaHCO3(aq), NaCO3 (aq) (see 01Table 1).

Groundwater temperatures range from 15 to 25 °C. The results of some ancillary nonisothermal runs indicate that changes in temperature within this range have a negligible influence on aqueous complexation and mineral dissolution-precipitation reactions. Therefore, calculations are isothermal at a temperature of 25 °C.

Initial and recharge (boundary) component concentrations are taken from Appelo (1994) and are listed in 02Table 2. The initial water (105 yr) is assumed to be brackish water as a result of the mixing of seawater with fresh water during the deposition of the overlying Marlboro clay, a brackish-water clay. The model assumes chemical equilibrium with respect to calcite. In the initial solution, the concentration of dissolved Ca2+ is such that the solution is in equilibrium with calcite. For a given pH, possible changes in the initial concentration of bicarbonate are compensated by changes in dissolved Ca2+ concentrations.

Prior information from Appelo (1994) is also used for initial guesses of model parameters and is incorporated into a least-squares objective function using the method described by Dai and Samper (2004). In the inversion, the computed concentrations are fitted to measured Na+, K+, Ca2+, Mg2+, HCO3, and pH spatial data.


Here we apply the inverse methodology for automatic identification of geochemical processes and estimation of model parameters. However, as pointed out by Carrera and Neuman (1986), whereas the formal term is “automatic parameter estimation,” not everything in the approach is truly automatic, nor should it be. Trial and error is carried out throughout the estimation procedure. The ultimate decision to select a good conceptual model and fit model parameters belongs to the modeler. We have used “preliminary” or “exploratory” trial-and-error inverse runs to estimate or fix some parameters before the formal inverse run. This procedure helps us to avoid overparameterization and reduces the risks of simultaneously estimating a large number of parameters.

The inverse problem of flow and reactive transport is solved in eight steps. 03Tables 3 and 044 list the results in each step, including the processes considered, parameters estimated, and the objective function values, as well as the contributions of each chemical species.

In step 1, the hydraulic head at the downstream boundary is fixed at 1.5 m based on Chapelle and Drummond (1983) and Appelo (1994). A flux is prescribed at the upstream boundary at 3.15 × 103 L/d, which was obtained from some preliminary inverse runs using head data. Appelo (1994) assumed that all the recharge is lost from the downstream half of the aquifer, across the leaky aquitard, to the upper aquifer, and there is no outflow at the downstream end of the aquifer. In order to test the plausibility of this assumption, a ground-water nodal leakage rate, Q, was estimated. In addition to Q, the transmissivities, T1 and T2, and the dispersivities (of the two parameter zones), α1 and α2, are estimated in this step. In this step only aqueous complexation and cation exchange reactions are considered. The result in this step showed that there is sensitivity to the outflow at the downstream end, and so we allow it to be nonzero and estimate its value in the steps that follow.

In step 2, calcite dissolution/precipitation at local chemical equilibrium is added to the geochemical system. In this step, a total of nine parameters are estimated: the two transmissivities, two dispersivities, nodal leakage rate, three selectivity coefficients, and CEC. The objective function reduces from 5194 to 4012, primarily due to the improvement on the fit of Na+ and pH, although Ca2+ and Mg2+ do not fit as well.

Step 3 is intended to evaluate the relevance of the processes postulated by Appelo (1994) and Plummer et al. (1994) to explain the relatively low observed pH values. This step is carried out in two stages. Step 3a considers proton surface complexation, while step 3b accounts for redox processes.

Accounting for proton surface complexation in step 3a requires adding XOH as a primary sorption species and considering the following surface complexation reactions:  
where log10K values for these reactions were obtained from Xu (1996). Since there are no measured data for the total density of sorption sites in Aquia aquifer materials, we took the value reported by Xu (1996) as prior information. Some exploratory inverse runs provided an estimate of the initial concentration of XOH of 5 × 10−2M. This is the total concentration of sorption sites, which is equal to the sum of the concentrations of surface sites XOH, XO, and XOH2+. We note that surface sites are not aqueous species, and, therefore, their concentrations are not comparable to those of dissolved species. Adsorbent specific surface area is 1000 m2/kg. A constant capacitance model is used for the surface electric potential (Appelo and Postma, 1993).
Step 3b considers a second mechanism for the release of protons. It is related to a source of CO2 coming from redox processes. Potential electron donors include organic matter of different redox states and methane. There are uncertainties on possible electron acceptors. The most obvious possibilities include oxygen, Fe(III), sulfate, and organic carbon (methano-genesis). Given the uncertainties in possible electron acceptors, our model does not consider any specific electron acceptor. Although iron and sulfate reduction could explain the observed low pH values, the lack of data on dissolved and mineral iron and sulfur phases prevents us from testing the relevance of these processes. Similar to what has been observed in other coastal aquifers (Dai and Samper, 2006) and following Appelo (1994), we consider that the CO2 source comes from methane oxidation. We model this process according to:  
where e denotes a free electron (Xu et al., 1999). In our numerical approach, redox reactions are modeled in terms of O2 (aq), which is treated as an aqueous primary species. Therefore, methane oxidation is modeled according to:  

Initial and boundary redox conditions were obtained after a few exploratory inverse runs, which led to a pE [pE = –log10(e)] of recharge water of 0 and an initial pE of −5. Such an initial pE corresponds to an almost complete reduction of sulfate in the system, which is consistent with the very low sulfate concentrations measured in the Aquia aquifer. In this step, the concentration of HCO3– in recharge water is estimated.

The objective function for step 3a (2660) is much smaller than that of step 3b (4536), as one can see in 04Table 4. The results of the inverse runs of steps 3a and 3b indicate clearly that proton surface complexation in step 3a is more important in explaining the measured data. Therefore, it can be concluded that, at least for the conditions considered in this step, the process of proton surface complexation is more plausible than methane oxidation. The possibility of both processes occurring simultaneously is explored in step 7.

Step 4 addresses the relevance of Mg-calcite dissolution. This mineral phase is added to the geochemical system considered in step 3a. For the time scales considered in the model, which cover a period of 105 yr, dissolution and precipitation of pure calcite and Mg-calcite can be generally assumed to take place at local chemical equilibrium. However, in this step, the possibility of reversible kinetically controlled dissolution precipitation of calcite and Mg-calcite is tested. These processes are modeled using kinetic rate laws with a specific surface of 20 m2/m3 for both minerals and rate constants of 1.36 × 10−12 mol m2 s−1 for calcite and 1.78 × 10−11 mol m2 s−1 for Mg-calcite. Modeling kinetic mineral dissolution and precipitation does not cause numerical problems, because initial time steps are small enough to ensure stability. Once equilibrium conditions are reached, time steps can be larger.

Two different volume fractions (fx) are considered: 0.1 in step 4a and 0.18 in step 4b. Neither of the two stages of step 4 lead to an improvement in the overall objective function compared to the results of step 3a (see 04Table 4). By accounting for the dissolution of Mg-calcite, the fit to most chemical species gets worse except for pH and HCO3. The results in 04Table 4 indicate that a Mg-calcite fraction of 0.18 leads to results that are slightly better than those obtained with a fraction of 0.1. The objective function decreases from 2864 for fx = 0.1 to 2785 for fx = 0.18. Surprisingly, the fit of Na+ is more sensitive to changes in the magnesium fraction of Mg-calcite than are the fits of Ca2+ and Mg2+. Since the results obtained in this step (4a and 4b) are worse than those in step 3a, Mg-calcite is eliminated in the following steps of the inverse problem, and calcite dissolution and precipitation at local chemical equilibrium is assumed henceforth.

Step 5 addresses the uncertainties in the initial concentrations of the chemical components in the aquifer. Initial concentrations of Na+, K+, Mg2+, and XOH are estimated by using the same model as step 3a. The objective function is reduced greatly from 2660 in step 3a to 1989 in step 5, mostly due to the improvement in the fit of Na+, Ca2+, and pH (see 04Table 4).

Step 6 follows a similar approach to step 5 and addresses the uncertainties in the concentrations of recharge waters. Boundary concentrations of Na+, K+, Ca2+, and HCO3 are estimated by using the same model as step 3a. The objective function shows a slight reduction from 2660 in step 3a to 2484 in step 6. Clearly, computed concentrations are less sensitive to recharge concentrations than they are to initial concentrations, possibly because prior estimates of recharge concentrations are close to their optimum values.

Step 7 aims at combining the conclusions of previous steps. It incorporates all the relevant chemical processes, including: cation exchange, proton surface complexation, methane oxidation, and calcite precipitation and dissolution. The following 17 parameters are estimated: 2 transmissivities, 2 dispersivities, groundwater discharge, CEC, 3 selectivities, initial concentrations of Na+, K+, Mg2+, and XOH, and recharge concentrations of Na+, K+, Ca2+, and HCO3. The objective function improves, decreasing to a value of 1661 04(Table 4). In this step, each chemical species contributes a relatively similar amount to the overall objective function (magnesium a little more), meaning that the fit is equitably good for all of them (see Fig. 3).

If initial and boundary waters are electrically neutral, computed aqueous solutions will remain neutral, because reactive transport models usually preserve charge. When either boundary or initial concentrations are estimated, there is no way to guarantee electroneutrality unless an additional constraint is added to the inverse model. Although the current version of our inverse code does not ensure electroneutrality, optimum estimates of initial and boundary concentrations in Aquia aquifer do not deviate from electroneutrality. Electrical imbalance calculated with EQ3/6 (Wolery, 1979) for recharge and initial concentrations is 2.03 × 10−12 of the total charge for the initial water and 5.83 × 10−13 for the recharge water. Therefore, both estimated initial and recharge concentrations reasonably satisfy the condition of electrical neutrality.

In step 8, the relevance of methane oxidation is revisited. By leaving out this reaction, the objective function increases from 1661 in step 7 to 1839 in step 8. This indicates that methane oxidation is a relevant process.

The best results of the inverse solution are therefore obtained in step 7. For the most part, estimated parameters, such as CEC, selectivities, and initial and recharge concentrations, are close to the values reported by Appelo (1994) (see 05Table 5). The relative differences between estimated results and those of prior information range from 3.8% to 90%. Estimated transmissivities are equal to 54 and 82 m2/d, which are within the range of measured values. Estimated dispersivities are also similar to those reported by Appelo (1994). Estimated CEC is 2.102 meq/100 g solid, which is smaller than the prior information (3.235 meq/100 g solid). This difference may be related to the different flow conditions used by Appelo for the downstream part, which will be discussed later.

Contrary to the trial-and-error method for parameter estimation, this inverse methodology provides objective measures for the reliability of parameter estimates, such as estimation variances, covariance and correlation matrices, and approximate confidence intervals. The correlation matrix of parameter estimates shows similar patterns to those of the column experiment addressed by Dai and Samper (2004). All three selectivity coefficients are positively correlated with each other and negatively correlated with respect to CEC. The correlation coefficients of CEC to KNa/K, KNa/Ca, and KNa/Mg are −0.12, −0.42, and −0.93, respectively. However, as stated in Dai and Samper (2004), this correlation does not influence the quality of estimated parameters when parameter prior information is considered in the estimation process. Uncertainties in estimated parameters, as measured by their standard deviations, are small (see 05Table 5).

The stepwise inversion method is useful for testing the plausibility of possible geochemical processes. The objective function is more sensitive to the different geochemical processes than to their parameters, which confirms that bias and uncertainty from inadequate conceptual models are much larger than those introduced from an inadequate choice of model parameter values (Hill, 1998; Neuman, 2003; Neuman and Wierenga, 2003; Meyer et al., 2004).

The best fit for different components after 105 yr is shown in Figure 3. Computed values match measured data for all six components with two exceptions. At the recharge area, the observed Ca2+ concentrations have a lot of scatter. The computed values go through the center of these scattered points. Computed Mg2+ concentrations are slightly larger than measured data at the downstream end of the flow path. This discrepancy could be due to having too little groundwater discharge through the downstream end. Our estimate of the leakage rate is equal to 96% of the total recharge. The remaining 4% discharges through the downstream boundary. Consequently, the water velocity in the downstream part of the aquifer is low, which leads to weak hydraulic dispersion and higher Mg2+ concentrations. Appelo (1994) used a noflow boundary at the downstream end and got an even greater discrepancy for concentrations of K+ and Mg2+.

Data on exchanged sodium were not used for the solution of the inverse problem and yet the process is represented very well in the results. The concentration of exchanged Na+ decreases with time as Ca2+ released by calcite dissolution replaces Na+. As a validation of the results of the inverse model, Figure 4 shows how the modeled exchange of Na+ compares to the measured data reported by Chapelle and Drummond (1983). This figure indicates that Ca2+ released by calcite dissolution replaces sodium within the exchange complex.

Figure 5 shows the time evolution of the cumulative amount of calcite precipitation (positive) or dissolution (negative). These results show that calcite is mainly dissolved in the Aquia aquifer. Downgradient from the recharge area, the increase in salinity follows with an increase in pH and HCO3 concentration and less dissolution of calcite.


Model parsimony and the identification of conceptual models can be assessed using identification criteria, as demonstrated by Carrera and Neuman (1986) in identifying flow models and by Samper and Neuman (1989) in inferring covariance functions for groundwater-quality data. These identification criteria are based on information theory and were developed within the context of maximum likelihood theory. Most of them minimize some measure of the closeness between the true (yet unknown) model and the proposed model. The first criterion is Akaike information criteria (AIC). The best model is that which minimizes the AIC criterion, defined as  
where S is the negative log-likelihood and Np the number of estimated parameters. This criterion lacks consistency, because the probability of choosing the wrong model does not go to zero as the sample size increases. The so-called modified Akaike information criterion (MAIC) is defined as  
where M is the number of data values used for solving the inverse problem. MAIC is consistent and therefore preferred. Hannan presented another modification of the AIC, the Hannan information criteria (HIC), which is defined as  
Our inverse methodology uses a dimensionless generalized least-squares objective function (see Dai and Samper, 2004), E(p), which, except for a constant, coincides with the negative log-likelihood function S (see Carrera and Neuman, 1986; McLaughlin and Townley, 1996).

Values of AIC, MAIC, and HIC have been calculated for all the alternative models considered for the Aquia aquifer. The number of data values (M) used for the inverse problem includes 28 observations for each one of the six chemical species (Ca2+, Mg2+, K+, Na+, HCO3, and pH), plus the number of data on prior information, which is equal to the number of parameters being estimated. Therefore, M = 6 × 28 + Np. One can see in 04Table 4 that the model of step 7 unequivocally minimizes the values of all three information criteria, thus indicating that this model is the best of all candidates.


Our stepwise inversion strategy consists of designing a list of plausible scenarios or steps, solving them in a stepwise manner, and selecting the one that provides the best solution to the inverse problem. This stepwise strategy allows us to test models of increasing complexity and avoid overparameterization. Identification of possible geochemical processes or scenarios is a subjective and hard-to-systemize task. One usually relies on previous works and available expertise.

There is always the possibility for other possible conceptual models. In formulating a given conceptual geochemical model, one should consider the availability of data for testing such a model. Our model of the Aquia aquifer did not consider iron reduction as a source of CO2, because no Fe data were available. As suggested by one of the anonymous reviewers, NETPATH-like mass-balance modeling could be used for the definition of the scenarios to be tested with inverse reactive transport modeling.

Our inverse modeling of the geochemical system in the Aquia aquifer focused on the main geochemical processes proposed in previous hydrochemical studies. There is, however, the possibility for other plausible conceptual models, such as those tested by Plummer et al. (1994). However, as recognized by Plummer et al. (1994), such models for redox processes can hardly be tested with an inverse model given the lack of relevant data on dissolved species and mineral phases. Identification of iron reduction processes cannot be achieved unambiguously unless dissolved and mineral iron data are available.

Our inverse modeling of hydrochemical data from the Aquia aquifer could be improved by extending the capabilities of our inverse code, INVERSE-CORE, to account for stable and radioactive isotopic data.

Chapelle and McMahon (1991b) found sulfate sources from confining beds to be relevant in the Black Creek aquifer in South Carolina. Their plausibility in the Aquia aquifer remains to be ascertained. However, in the Aquia aquifer, water leaks out of the aquifer, which makes it unlikely that sulfate leaches into the aquifer against water flow by backward diffusion.

The inverse geochemical model of the Aquia aquifer could be improved by reassessing the role of Mg mineral sources, such as low Mg-calcite, and testing the relevance of kinetically controlled decomposition of organic matter with a redox state near zero.

Another source of uncertainty in our model comes from the assumption of constant hydrologic conditions over the past 100 k.y. Noble gas temperatures show the presence of water that infiltrated under conditions that, during the Last Glacial Maximum, were 9 °C colder than at present. Noble gas temperatures correlate with chloride concentrations, suggesting that chloride variations in the aquifer constitute a climate signal (Aeschbach-Hertig et al., 2002). Stable isotope ratios, however, do not provide a clear record of past climatic changes in the Aquia aquifer (Aeschbach-Hertig et al., 2002). Given the lack of conclusive evidence on time-varying conditions for the downstream-end boundary, our estimated discharge rate through the downstream-end boundary must be considered as an average value. Discharge should be smaller (larger) under higher (lower) sea-level conditions.


We have presented a stepwise inverse modeling methodology that is capable of testing and comparing possible conceptual geochemical models or scenarios and selecting the best according to available hydrochemical data, which do not necessarily have to be located along a flow path. The inverse methodology has been applied to modeling the geochemical processes involved in the displacement of seawater by fresh water in the coastal Aquia aquifer during the last 105 yr. The inverse model accounts for a wide range of chemical processes, such as aqueous complexation, acid-base and redox reactions, cation exchange, proton surface complexation, and mineral dissolution and precipitation.

Identification criteria derived from information theory were used to select the best among ten candidate conceptual models. Results show that cation exchange and calcite dissolution are the dominant geochemical process. Although both proton surface complexation and methane oxidation can explain the relatively low measured pH values, the former is more effective in fitting measured hydrochemical data. The results of the inverse model indicate that both processes may have occurred simultaneously in Aquia groundwater. On the other hand, the inverse model rejects the scenario of Mg-calcite dissolution and precipitation.

We obtained optimum estimates of transmissivities, dispersivities, groundwater discharge, selectivity coefficients, CEC, and initial and recharge concentrations, which lead to a fit of measured concentrations of Ca2+, Mg2+, K+, Na+, HCO3, and pH that improves upon the results obtained by Appelo (1994). The relative differences between our estimates and prior estimates range from 3.8% to 90%. Uncertainties in initial concentrations are larger than uncertainties in concentrations of recharge water. Estimated values of initial concentrations indicate that the initial water in the Aquia aquifer was not pure seawater but a mixture of seawater and fresh water. In the final model, 96% of recharge water leaks through the confining layers, and a small but significant part of the flow (4%) makes it through the downstream end of the aquifer.

This work was supported by the Full-scale Engineered Barriers Experiment in crystalline host rock (FEBEX) Project, funded by Spanish Nuclear Waste Company (ENRESA) (contracts 70323 and 770045) and the European Union (Projects FI4W-CT95-0008 and FIKW-CT2000-0016) within the Nuclear Fission Safety Programme; the Spanish Ministry of Science and Technology (CICYT HID98-282); and the University of La Coruña through a research scholarship awarded to the first author. We are grateful to Tianfu Xu and Francis H. Chapelle for providing basic data as well as sharing their interesting thoughts about the Aquia aquifer.