Use of geochemical modeling to evaluate the hydraulic connection
of aquifers: a case study from Chianan Plain, Taiwan
Hsueh-Yu Lu&Tsung-Kwei Liu&Wen-Fu Chen& Tsung-Ren Peng&Chung-Ho Wang&Meng-Hsun Tsai& Tai-Sheng Liou
Abstract Aquifers are generally composed of highly permeable layers that can conduct a considerable amount of groundwater. Traditionally, aquifer units are correlated through the concept of lithostratigraphy. For low-permeable aquifers, it is difﬁcult to deﬁne the spatial distribution of hydrogeological units, and this study attempts to use geochemical modeling to identify the groundwater ﬂow paths in an area of Taiwan. Multiple geochemical analyses, including groundwater chemistry; stable isotopic composi-tions of hydrogen, oxygen and carbon; and radiocarbon contents were performed. Using these parameters as the constraints of geochemical models, the hydraulic
connec-tion was examined between pairs of possibly interlinked wells along four selected cross sections, and the conceptual groundwater model was accordingly established. The resultant model suggests that the hydraulic connection between aquifers should be correlated with the concept of chronological stratigraphy, especially for low-permeable, unconsolidated aquifers. Using Darcy’s law, the hydraulic conductivities of theﬁne-sand aquifers were estimated to be between 3.14 × 10−5 and 1.80 × 10−4 m/s, which are roughly one order of magnitude higher than those derived by in situ pumping tests. The substantial extraction of groundwater over a long period in the studied area could accelerate groundwaterﬂow, leading to an overestimation of the aquifer permeability.
Résumé En règle générale, les aquifères sont constitués de couches très perméables capables de conduire des quantités d’eau considérables. Les unités aquifères sont tradition-nellement corrélées d’un point de vue lithostratigraphique. Dans le cas des aquifères à faible perméabilité, il s’avère difﬁcile de déﬁnir une distribution spatiale des unités hydrogéologiques, et cette étude tente d’utiliser la modéli-sation géochimique pour identiﬁer les cheminements des eaux souterraines dans un secteur de Taiwan. Plusieurs analyses géochimiques ont été réalisées ; elles incluent la chimie des eaux souterraines, les isotopes stables de l’hydrogène, de l’oxygène et du carbone et l’abondance en carbone 14. En utilisant ces paramètres comme contraintes pour les modèles géochimiques, les liaisons hydrauliques entre les puits potentiellement interconnectés ont été étudiées deux par deux, le long de quatre coupes présélec-tionnées, et le modèle conceptuel a été établi en consé-quence. Ce modèle résultant suggère que les connexions hydrauliques inter-aquifères devraient être corrélées avec les concepts stratigraphiques, surtout pour les aquifères non consolidés à faible perméabilité. Les perméabilités dans les aquifères de sablesﬁns ont été estimées par la loi de Darcy entre 3.14×10−5 et 1.80×10−4 m/s, soit environ un ordre de grandeur au-dessus de celles issues des tests de pompage sur site. L’exploitation substantielle des eaux souterraines sur une longue période dans la zone d’étude peut accélérer les écoulements souterrains, menant à une surestimation de la perméabilité de l’aquifère.
Received: 21 September 2006 / Accepted: 9 July 2007 Published online: 24 August 2007
Electronic supplementary material The online version of this article (doi:10.1007/s10040-007-0209-6) contains supplementary material, which is available to authorized users.
H.-Y. Lu ())
Department of Earth and Environmental Sciences, National Chung Cheng University,
168, University Rd.,
Min-Hsiung, Chia-Yi, Taiwan e-mail: email@example.com Tel.: +886-5-2720411 Fax: +886-5-2720807 T.-K. Liu
Department of Geosciences, National Taiwan University, Taipei, Taiwan
Department of Tourism Management,
Chia Nan University of Pharmacy and Science, Chia-Yi, Taiwan
Department of Soil and Environmental Sciences, National ChungHsing University,
Taichung, Taiwan C.-H. Wang
Institute of Earth Sciences, Academia Sinica, Taipei, Taiwan
Resumen Los acuíferos están generalmente compuestos de capas muy permeables que pueden conducir una cantidad considerable de agua subterránea. Tradicional-mente, se ponen en correlación las unidades acuíferas a través del concepto de litoestratigrafía. Para los acuíferos de baja-permeabilidad, es difícil de deﬁnir la distribución espacial de unidades hidrogeológicas y este estudio intenta usar modelamiento geoquímico para identiﬁcar las direcciones delﬂujo de agua subterránea en un área de Taiwán. Se realizó el análisis geoquímico múltiple, incluyendo la química del agua subterránea, la composi-ción de isótopos estables de hidrógeno, oxígeno y carbono, y el contenido de radiocarbono. Usando estos parámetros como limitantes de modelos geoquímicos, la conexión hidráulica se examinó entre pares de pozos posiblemente interrelacionados, a lo largo de cuatro cortes transversales seleccionados y de acuerdo con esto se estableció el modelo conceptual del agua subterránea. El modelo resultante sugiere que la conexión hidráulica entre los acuíferos deba interrelacionarse con el concepto de estratigrafía cronológica, sobre todo para los acuíferos sin consolidar de baja permeabilidad. Usando la ley de Darcy, se estimaron las conductividades hidráulicas de los acuíferos de arena ﬁna obteniendo valores entre 3.14× 10−5y 1.80×10−4m/s, las cuales son aproximadamente un orden de magnitud mayor, que aquéllos obtenidos in situ por las pruebas de bombeo. La extracción considerable de agua subterránea durante un periodo largo en el área estudiada, podría acelerar el ﬂujo de agua subterránea, llevando a una sobrestimación de la permeabilidad en los acuíferos. Keywords Aquifer correlation . Chronostratigraphy . Geochemical modeling . Groundwater age . Taiwan
Determining the spatial distribution of hydrogeological units is at the core of establishing a conceptual model of a regional groundwaterﬂow system. A number of analytical methods can be used to assess the lateral and vertical connections of hydrogeological units. In general, aquifers are composed of lithologic units of high hydraulic conduc-tivity, such as porous basaltic layers, fractured sandstone, and unconsolidated sand–gravel formations (e.g., Sugarman and Miller1997; Hiatt et al.2003). These aquifers can be correlated by means of their lithology, sedimentary structures, and thickness. To integrate surface geophysical and cross-hole geophysical surveys, high-resolution hydro-geological information can be determined in remote and noninvasive ways (e.g., Shtivelman and Goldman 2000; Lubczynski and Roy 2003). However, these methods are all conducted on the basis of lithostratigraphic relation-ships. For poorly conductive aquifers and highly deformed layers, aquifer connectivity is rarely consistent with stratigraphic correlation, and hydraulic potential may vary greatly both in space and time to hinder the determination of aquifer layering. Under these circumstances, a suitable technique is essential to genuinely track the groundwater
ﬂow path. In this study, radiocarbon (14C) groundwater
dating is utilized as a natural tracer.
The application of 14C dating to groundwater systems encounters two major problems: (1) determining the initial concentration of radiocarbon for the recharge water, and (2) quantifying the mass transfer of radiocarbon between soil and groundwater due to geochemical reactions, including evaporation and groundwater mixing (Plummer et al.1994). Speciﬁcally, the measurement of 14C concen-tration should take account of geochemical reactions affecting 14C activity in groundwater systems, yielding a “corrected” 14C age. This is generally approached by
inverse geochemical modeling between two points within a water body along the ﬂow path(s), if sufﬁcient data are available, including chemical composition of groundwater, isotopic composition of groundwater, and mineralogy of sediments (Buckau et al.2000; Van der Kemp et al.2000). The modeled (corrected)14C age of groundwater has been widely used to establish regional groundwaterﬂow systems (e.g., André et al.2005; Hidalgo and Cruz-Sanjulián2001) and groundwater recharge patterns (e.g. Bajjali 2006; Cronin et al.2005; Zuber et al.2000).
Many programs have been developed to interpret the isotopic and geochemical evolution of groundwater with inverse geochemical modeling. The geochemical mass-balance and mixing computer code NETPATH was employed in this study. NETPATH is an interactive program developed by the US Geological Survey (USGS) for inversely modeling net chemical and isotope mass-balance reactions between initial andﬁnal waters along a hydrogeological ﬂow path (Plummer et al. 1994). NET-PATH computes the carbon isotopes (δ13C) ofﬁnal water as a Rayleigh distillation problem using the equations of Wigley et al. (1978) for the modeled mass transfer. NETPATH is widely used to determine residence time of groundwater based on carbon isotopes.
The Taiwan Government initiated the nationwide Taiwan Groundwater Monitoring Network (TGMN) in 1992 to investigate Taiwan’s groundwater resources. So far, this project has installed more than 900 monitoring wells and provides versatile hydrogeological information for exploi-tation and conservation of groundwater resources. The study area is in the North Chianan Plain Groundwater District (NCPGD), which is one of the regional monitor-ing networks of the TGMN. NCPGD is located along the western coast of southern Taiwan, which is a foreland basin bounded by the Western Foothills on the east (Fig. 1). There are six major river systems in the study area. From north to south they are: Bajhang Creek, Jishuei Creek, Jiangjyun Creek, Cigu Creek, Cengwun Creek, and Yanshuei Creek, which are geographically classiﬁed as meandering due to gentle regional slopes of tens to a hundred centimeters per kilometer and low bed loads. The Western Foothills began uplifting during the Pleistocene when the collision of the Eurasian Plate with the Philippine
Sea Plate occurred, and the basin now correspondingly subsides at a rate of 5–11 mm/year, increasing from north to south and from east to west. The sediments, accumulated over the Holocene period, are 80-m thick in the north of the NCPGD and increase to >100 m in the coastal area to the southern end of the NCPGD. The total thickness is more than 300 m. In general, the aquifers that developed in the basin were deposited in environments of Holocene– Pleistocene continental shelf to littoral sea and are mainly composed of silt and clay sediments. These ﬁne-grained aquifers are caused by the dominance of Pliocene and Miocene mudstone in the Western Foothills, which supplied most of the sediments for the NCPGD (Fig.1). Therefore, the NCPGD is generally short of groundwater resources, as the spatial and temporal distributions of hydraulic potential are highly variable. It is difﬁcult to determine the extent of aquifers by conventional methods owing to the low density of monitoring wells.
The aquifers in the NCPGD were recharged mainly along the western front of the Western Foothills, which is composed of semiconsolidated Pleistocene sandstone with minor conglomerate deposited on the river terrace (Fig.1). The top aquifer is also naturally replenished by surface
water from precipitation and streams. As a result, groundwater preferentially discharges from the Western Foothills on the east side to the western coastal area.
A high concentration of arsenic in the groundwater is a notorious phenomenon in the NCPGD (Lu 1990). From the beginning of the twentieth century, citizens extracted shallow groundwater from the top aquifer deposited during the Holocene high sea-level period. The water was of poor quality for daily usage but was low in arsenic. Since the 1950s, a shortage of water resources has encouraged people to install wells of 100–200 m deep into strata deposited in the last glacial period (18,000 years BP). In contrast to the waters of the top aquifer, the water from the lower aquifer is less salty but has a higher arsenic concentration of 0.1–2.5 mg/, as shown in water-quality monitoring data from the Inquiry and Supply System of Water Resources Database (Water Resources Agency2006). The spatial distribution of arsenic content is apparently different to that in Bangladesh (Van Geen et al. 2003), where arsenic concentrations are higher in aquifers closer to the surface. The problems associated with long-term consumption of the arsenic-tainted water are evidenced in the population, with various severe diseases Fig. 1 Geological map of North Chianan Plain Groundwater District (NCPGD) showing the 22 sampled sites (N01-N22) and four sections (I–IV) along the groundwater ﬂow path used for geochemical modeling
including skin problems such as “black-foot disease.” However, the problem is no longer epidemic, and only a few new patients have been reported in recent years.
Groundwater samples were collected for geochemical analysis at different depths from 22 locations during 2004. The sample numbers are preﬁxed by a letter N, corresponding to NCPGD; the remainder consists of two digits for site number followed by one digit for screen number, such as N01-1 and N02-2. The smaller the screen number, the shallower the sample. Prior to collecting samples, at least three volumes of water in the well casings were extracted at a reasonably low rate and drained away from the well area to prevent contamination from stagnant water. During well purging, water quality was frequently recorded. The samples were collected when dissolved oxygen (DO), temperature, pH, electrical conductivity (EC), Eh (redox potential), and turbidity remained stable. Once suitable samples were obtained, they were ﬁltered through 0.45-μm glass ﬁber papers and acidiﬁed with ultrapure nitric acid to pH 2 for subsequent laboratory analysis. Cations were determined by inductively coupled plasma–optical emission spectroscopy (ICP-OES), whereas anions were analyzed by ion chromatography.
For analysis of stable isotopic composition (oxygen, hydrogen, and carbon) and radiocarbon activity, groundwa-ter samples were collected from newly constructed monitor-ing wells just after well completion and performance tests during the period 2000–2002. The reason for collecting samples in this manner was to diminish the possibility of isotopic alteration by chemical and biological reactions after well construction. The oxygen isotopic ratios were deter-mined using isotopically equilibrated CO2 gas generated using CO2-H2O equilibration (Epstein and Mayeda1953), whereas the hydrogen isotopic ratios were analyzed using hydrogen gas produced by reducing water with zinc (Coleman et al.1982). Both analyses were conducted with isotope ratio mass spectrometers at the Isotope Hydrology Laboratory of Academic Sinica, Taiwan. The oxygen and hydrogen isotopic ratio results are reported usingδ-notation as per mil (‰) relative to an international Vienna Standard Mean Ocean Water (VSMOW) standard and normalized on scales such that the δ18O and δD of Standard Light Antarctic Precipitation (SLAP) are −55.5 ‰ and −428 ‰, respectively. In this study, the analytical errors (1σ) are 0.1 ‰ for δ18
O and 1.5‰ for δD.
For carbon isotope analysis, dissolved inorganic carbon (DIC) was precipitated from water, in theﬁeld, by adding BaCl2solution. The samples were then transported to the
Environmental Isotope Laboratory, National Taiwan Uni-versity, Taiwan, for measurement. The carbon in water samples was extracted and synthesized into benzene (92% carbon). Radiocarbon activity, given in percent modern carbon (pmc), was subsequently determined by a liquid scintillation counting technique. Carbon-13 (13C) was analyzed by mass spectrometry, and the resultant data
are expressed in δ-notation relative to Vienna Pee Dee Belemnite (VPDB) standard. The analytical errors are 0.4–0.7 pmc for 14C content and 0.1 ‰ for δ13C, respectively. Moreover, shell fragments and carbon sam-ples were also collected from drilling cores. The samsam-ples were sent to Beta Analytic Inc., USA, and Waikato Radiocarbon Dating Laboratory, New Zealand, for 14C age dating. The results can be retrieved from the Hydro-geology Database (Central geological Survey2006).
The results of chemical analysis for 68 groundwater samples from 22 sites (N01-N22) in the NCPGD are listed in electronic supplementary material (ESM) Table1. Generally speaking, the groundwaters in the NCPGD are typically high in EC, ranging from 537 to 103,600μS/cm, and the Piper diagram shows that most of the samples are distributed below the seawater mixing line (Fig. 2). The samples collected from the top aquifer near the recharge area, at sites N03, N09, and N22, are generally identiﬁed as Ca-HCO3hydrochemical facies resulting from surface-acidic water dissolving carbonate. The groundwaters (sites N02, N08, N16, and N21) evolved to Na-HCO3 facies when ﬂowing down gradient. This evolutionary trend suggests a geochemical process resulting from cation exchange (Appelo1994; Güler and Thyne 2004):
2þþ Na X ! 1
2Ca X2þ Na
where X represents exchanging substrates. This reaction explains the observed increase in Na+ concentration without an associated increase in Cl− concentration along theﬂow path. Based on the above reaction, the aquifer is primarily dominated by Na-rich clay. It represents a process whereby a brackish aquifer isﬂushed with fresh carbonate water. As aforementioned, the aquifers in the NCPGD are mostly composed of poorly permeable materials and permit a very low recharge rate. The syndepositional connate seawater trapped within aquifers was only partially expelled in the coastal area, as shown by evidence of seawater in samples from sites N01, N02, N04, N06, N07, N13, N14, and N15. These samples show transitional geochemical characteristics between Na-HCO3 facies and Na-Cl facies
in the Piper diagram. In addition, Ca-Mg-HCO3-SO4facies
groundwaters are found in samples N01-1, N02-1, and N22-1, which may involve gypsum dissolution owing to the high Ca and high sulfate but low sulﬁde hydrochemical facies. The high sea level of the mid-Holocene submerged most of the area of the NCPGD; the consequent develop-ment of lagoon sedidevelop-mentation provided an ideal environ-ment for the precipitation of gypsum.
Scatter plots of the major ion components of the groundwaters also demonstrate the two stages of geo-chemical evolution, as illustrated in the Piper diagram, which are controlled by cation exchange between ﬂuid and clay mineral and seawater freshening (or freshwater ﬂushing). Compared with a nearly linear relationship
between Na+ and Mg2+ vs. Cl− (Fig. 3a and 3b, respectively), the plot of Ca2+ vs. Cl− shows a scatter cluster where Cl−concentration is under 100 mmol/l and a linear relationship on the opposite side (Fig. 3c). It demonstrates that Ca-carbonate dissolution is the domi-nant reaction in a freshwater system, whereas that of Mg-carbonate is not signiﬁcant. Subsequently, the resultant groundwater ﬂushes the brackish aquifer to show geo-chemical characteristics of seawater freshening.
Hydrogen and oxygen stable isotopes
In this study, hydrogen and oxygen stable isotopes were utilized as a constraint to calculate the mixing proportion of seawater. As mentioned, seawater plays an important role in the geochemical composition of the NCPGD. Having a better estimation of the mixing proportion of seawater is essential to simulate geochemical evolution. However, only six samples from two sites (N03 and N19) from the NCPGD were analyzed for hydrogen and oxygen stable isotopic composition (ESM Table 1). Fortunately, there are 27 analyses available for stable isotopes and basic geochemical properties from the South Chianan Plain Groundwater District (SCPGD) located adjacent to
the NCPGD on the southern boundary (Yanshuei Creek) (Chen et al. 2004). The correlation between the mixing proportion of seawater and geochemical properties for the SCPGD can be established. By applying this relationship to the NCPGD, the mixing proportion of seawater is derived for the subsequent geochemical modeling.
Relationship between isotopes and mixing proportionThe stable isotopic compositions of water from the SCPGD apparently show a linear relationship between two end members: seawater (δD≈0 ‰, δ18O≈0 ‰) and freshwater (δD≈−42 ‰, δ18O≈−7 ‰) (Fig.4). On the lighter isotope side, the scatter cluster represents isotopic variation of fresh groundwaters. The linear relationship shows a process of seawater mixing and is utilized to estimate the mixing proportion of seawater. To exclude the data having hydrogen isotopes lighter than−42.0 ‰ (pure freshwater), the isotope data from the SCPGD were ﬁrst regressed to pass through the point ofδD=0 ‰ and δ18O=0‰, which is assumed to be the representative isotopic composition of pure seawater. The resulting equation, with correlation coefﬁcient (r2) of 0.972, can be written as (Fig.4):
%D ¼ 5:57 %18O ð2Þ
Table 1 Isotopic composition of groundwaters from the South Chianan Plain Groundwater District (SCPGD) (Peng2001) Sample
Depth (m) δD (‰) δ18O (‰) Electrical conductivity
(EC) (μS/cm) ProjectedδD (‰)b Unadjusted mixing proportion of seawater (%)c Adjusted mixing proportion of seawater (%)d 1-1 53 −26.9 −4.5 19,720 −26.8 35.9 24.1 1-2 230 −40.2 −6.6 2,660 −40.1 4.3 0.0 2-1 1 −39.0 −6.6 857 −38.9 7.1 0.0 3-1 29 −48.5 −7.2 2,150 −42.0 0.0 0.0 3-2 203 −30.5 −5.3 3,700 −29.5 27.3 13.9 4-1 25 −39.6 −6.3 2,062 −38.4 5.8 0.0 4-2 106 −46.0 −7.1 554 −42.0 0.0 0.0 5-1 44 −3.7 −0.5 77,200 −3.7 91.2 89.6 5-2 143 −2.4 −0.4 48,900 −2.4 94.3 93.2 6-1 7 −49.1 −8.3 1,648 −42.0 0.0 0.0 7-1 57 −2.9 −0.4 65,300 −2.9 93.1 91.8 7-2 155 3.0 0.4 41,800 0.0 100.0 100.0 8-1 13 −4.0 −1.8 8,810 −4.2 90.0 88.2 8-2 108 −16.5 −3.1 29,800 −16.5 60.6 53.3 8-3 199 −3.6 −1.5 39,300 −3.8 91.0 89.4 9-1 6 −35.3 −6.8 885 −35.4 15.5 0.0 9-2 143 −41.0 −7.1 1,159 −40.9 2.3 0.0 10-1 20 1.0 0.1 66,900 0.0 100.0 100.0 11-1 2 −33.9 −6.3 1,839 −33.9 19.0 4.1 12-1 6 −45.4 −7.1 9,760 −42.0 0.0 0.0 13-1 13 −29.5 −5.4 3,410 −29.5 29.6 16.6 13-2 158 −33.8 −6.0 7,890 −33.8 19.4 4.5 13-3 199 −28.9 −5.6 9,410 −29.0 30.9 18.2 14-1 4 −42.0 −7.0 1,425 −42.0 0.0 0.0 15-1 19 −11.4 −2.1 41,200 −11.4 72.8 67.8 15-2 73 −23.4 −4.3 27,800 −23.4 44.1 33.8 15-3 212 −19.6 −3.7 30,100 −19.6 53.1 44.5 a
The sample locations are shown in electronic supplementary material (ESM) Fig.1. The sample numbers consist of a number for site number followed by one digit for screen number
The projectedδD values were calculated by projecting the original δD values on the isotopic mixing line (Eq.2) c
The unadjusted mixing proportions of seawater were calculated based on the interpolation relationship between presumed pure seawater and pure freshwater (Eq.3)
The original data points were perpendicularly projected on the regression line to recalculate the projected isotopic compositions (Table 1). As shown in Fig. 4, the isotopic compositions of δD=−42.00 ‰ and δD=0 ‰ were selected as the representative values of pure freshwater and pure seawater, respectively. The values beyond this range were cut off in the following calculations. Assuming a linear relationship of stable isotopic composition be-tween pure seawater and pure freshwater, the mixing proportion of seawater can be interpolated as:
UMP ¼ 100 ð%D þ 42:0Þ
where UMP is the unadjusted mixing proportion of seawater in Table 1. At this point, the calculated mixing proportion is based on a roughly selected isotopic composition of freshwater. It needs to be further adjusted.
Relationship between EC and mixing proportionTo compare with other basic geochemical properties, the EC value is found to have a better correlation with the unadjusted mixing proportion of seawater (Fig. 5). How-ever, Fig.5shows three distinct groups. The EC values on the left side of theﬁgure almost remain constant and keep low, whereas unadjusted proportions increase. It means that this group can be considered as isotopic variation of Fig. 2 Piper diagram of water chemistry of the North Chianan Plain Groundwater District (NCPGD) aquifers. Black dots represent the 68 groundwater samples taken in 2004. Electrical conductivity (EC) concentration is represented on a proportional scale relative to the sizes of the circles shown in the legend
pure freshwater. The group on the right of theﬁgure with high EC shows a scatter relationship with unadjusted mixing proportion of seawater. This may be caused by complicated chemical reactions in high-ionic-strength systems that signiﬁcantly modify geochemical properties, suggesting that geochemical modeling should not be applied in such a system. Data points in the middle of the ﬁgure show a linear positive correlation, and the lightest isotopic composition is reselected as the represen-tative isotopic composition of pure freshwater (δD=−35.4 ‰ and δ18
O=−6.4 ‰). The adjusted mixing proportion of seawater is accordingly rescaled, as listed in Table 1. Consequently, the adjusted mixing proportions were
correlated with EC values to obtain a regressed equation with a correlation coefﬁcient (r2) of 0.897:
EC ¼ 623:46 þ 600:58 AMP ð4Þ
where EC is electrical conductivity inμS/cm and AMP is the adjusted mixing proportion of seawater. This equation was applied to calculate the mixing proportion of seawater in the NCPGD, which was used as a constraint in the NETPATH model (Table 2).
Geochemical reaction system
In this study,13C and14C are regarded as major criteria for obtaining a plausible geochemical model using NET-PATH. Their analytical results for 67 samples from the NCPGD are listed in ESM Table 1. To establish a geochemical reaction system concentrating on the mass balance of carbon isotopes is the ﬁrst step for inverse geochemical modeling. As indicted in the previous discussions, the groundwater ﬂow reaction system in the NCPGD includes three major geochemical processes: carbonate dissolution/precipitation, organic and inorganic carbon dissolution/degassing, and sulfur-related redox reactions (André et al. 2005). The former two processes are commonly dominant in most groundwater systems, and the last one is important in brackish aquifers with a high level of sulfur, which may be ignored in a freshwater system. Accordingly, the elements used in NETPATH to constrain the mass-balance solution include C, S, Ca, Mg, Na, Fe, and redox state (Fig. 6). The redox state is a constraint that must be included in the redox reactions, as suggested by Plummer et al. (1994). The model also employed δD or δ18O, as listed in Table1, to determine the mixing proportion of seawater. Based on the selected constraints, the reactant phases, including calcite, dolo-mite, cation exchange between Ca2+ and Na-X, CH2O,
CH4, CO2, gypsum, siderite, geothite, pyrite, and H2S,
were used in NETPATH (Fig.6).
Presumed chemistry of seawater and infiltration
For modeling top aquifers, surface inﬁltration should be considered important. However, the geochemical proper-ties of inﬁltration water are difﬁcult to measure; surface inﬁltration is calculated by using PHREEQC. PHREEQC is a USGS computer code for aqueous geochemical calculations that can solve speciation, inverse, reaction-path, and 1-D advective reaction-transport modeling problems (Parkhurst and Appelo 1999). Assuming pure water equilibrated with sedimentary calcite, and assuming 10−3.5 atm soil CO2 gas in an open system at 25°C, the
geochemical properties have been simulated, as listed in Table3. According to the model of Ingerson and Pearson (1964), the carbon isotopes in inﬁltration water are contributed nearly equally from sedimentary carbonate and soil CO2 gas and can be roughly estimated by the Fig. 3 Chemical composition of major elements vs. Cl−in North
Chianan Plain Groundwater District (NCPGD): Na+(a) , Mg2+(b) , and Ca2+(c) . The best-ﬁt lines and their equations are shown
mean value of these two sources. The default values of carbon isotopes in NETPATH soil CO2gas are traditionally designated as −25 ‰ for δ13C and 100 pmc for 14C, whereas marine carbonate is assumed to have a δ13C content of 0 ‰ and14C content of 0 pmc. However, there are many 14C ages available for subsurface sediment
samples from variable depths in the NCPGD (Hydro-geology Database, Central Geological Survey 2006). The top aquifers at depths of about 50 m were generally deposited in the last 6,000 years, which is roughly equivalent to14C content of 50 pmc instead of the default value of 0 pmc. The carbon isotopic compositions of in
ﬁl-Fig. 5 Correlation of electrical conductivity (EC) and unadjusted mixing proportion of seawater. The group on the lower left end of the regression line is considered as isotopic variation of pure freshwater. The scattered group on the right-hand side with high EC suggests that geochemical modeling may not be applicable in high EC waters. AMP adjusted mixing proportion of seawater
Fig. 4 δD–δ18O plot of groundwaters from South Chianan Plain Groundwater District (SCPGD) aquifers. The regression line is shown as the solid line through presumed isotopic compositions of pure seawater and pure freshwater (marked with open squares)
tration water are subsequently estimated as 75 pmc for14C and−12.5 ‰ for δ13C, on average (Table3). Furthermore, the geochemical properties and carbon isotopes of seawater are also essential for the mixing model of a brackish
aquifer. In this study, seawater composition proposed by Goldberg et al. (1971) was employed (Table 3); the δ13C was traditionally designated as 0 ‰ and 14C content was assigned to sediment14C concentration.
Table 2 Derived isotopic compositions of groundwaters from the North Chianan Plain Groundwater District (NCPGD) as a constraint for NETPATH model
Well Electrical conductivity
(μS/cm) DerivedδD (‰) Derivedδ18O (‰) Well Electrical conductivity(μS/cm) DerivedδD (‰) Derivedδ18O (‰)
N01-1 900 −35.2 −6.3 N11-3 630 −35.4 −6.4 N01-2 1,880 −34.6 −6.2 N12-1 699 −35.3 −6.3 N01-3 3,670 −33.6 −6.0 N12-2 829 −35.3 −6.3 N01-4 729 −35.3 −6.3 N12-3 998 −35.2 −6.3 N02-1 665 −35.4 −6.3 N13-1 8,920 −30.5 −5.5 N02-2 1,221 −35.0 −6.3 N13-2 4,420 −33.1 −6.0 N03-1 516 −35.4 −6.4 N14-1 76,700 0.0 0.0 N03-2 539 −35.4 −6.4 N14-2 16,500 −26.0 −4.7 N03-3 660 −35.4 −6.3 N15-1 3,800 −33.5 −6.0 N04-1 61,100 0.0 0.0 N15-2 31,900 −17.0 −3.0 N04-2 1,340 −35.0 −6.3 N15-3 981 −35.2 −6.3 N04-3 878 −35.2 −6.3 N15-4 1,537 −34.8 −6.3 N05-1 1,001 −35.2 −6.3 N16-1 851 −35.3 −6.3 N05-2 458 −35.4 −6.4 N16-2 918 −35.2 −6.3 N05-3 514 −35.4 −6.4 N16-3 725 −35.3 −6.3 N05-4 507 −35.4 −6.4 N16-4 687 −35.3 −6.3 N06-1 4,650 −33.0 −5.9 N17-1 2,970 −34.0 −6.1 N07-1 50,500 −6.0 −1.1 N17-2 16,290 −26.2 −4.7 N07-2 9,210 −30.3 −5.4 N17-3 1,260 −35.0 −6.3 N07-3 4,540 −33.1 −5.9 N17-4 1,082 −35.1 −6.3 N07-4 1,011 −35.2 −6.3 N18-1 1,007 −35.2 −6.3 N08-1 469 −35.4 −6.4 N18-2 1,167 −35.1 −6.3 N08-2 3,830 −33.5 −6.0 N18-3 713 −35.3 −6.3 N08-3 1,096 −35.1 −6.3 N18-4 667 −35.4 −6.3 N08-4 3,220 −33.9 −6.1 N19-1 98,300 0.0 0.0 N09-1 1,008 −35.2 −6.3 N19-2 12,940 −28.1 −5.0 N09-2 725 −35.3 −6.3 N19-3 1,485 −34.9 −6.3 N09-3 513 −35.4 −6.4 N20-1 67,900 0.0 0.0 N10-1 53,400 −4.3 −0.8 N20-2 3,790 −33.5 −6.0 N10-2 41,200 −11.5 −2.1 N20-3 11,280 −29.1 −5.2 N10-3 14,730 −27.1 −4.9 N21-1 934 −35.2 −6.3 N10-4 932 −35.2 −6.3 N21-2 1,024 −35.1 −6.3 N11-1 1,094 −35.1 −6.3 N21-3 874 −35.2 −6.3 N11-2 17,600 −25.4 −4.6 N22-1 1,240 −35.0 −6.3
Fig. 6 Geochemical process-es governing the carbon isotopes of groundwaters from North Chianan Plain Ground-water District (NCPGD) (André et al.2005)
NETPATH models and flow system
Four sections (sections I–IV) along the groundwater ﬂow path were modeled in this study (Fig. 1). The hydraulic connection between pairs of possibly interlinked monitor-ing wells was examined by the NETPATH model. There are multiple solutions for the geochemical reactions, as shown in Fig.6, due to a greater number of reactant phases being employed to ﬁt a lesser number of constraints. Generally, the models with calculatedδ13C roughly ﬁtting the observed value were considered as the possible geochemical interpretation for hydraulically connected aquifers. The best-ﬁt model was further conﬁrmed with reasonable travel time and mineral mass transfer. The simulated results are summarized in Appendix. However, a plausible model may not be available. It suggests that two modeled aquifers are not located on the sameﬂow path. In addition, NETPATH uses WATEQ4F to model the thermo-dynamic speciation of major and important minor inorganic ions (Plummer et al. 1994). Although WATEQ4F shows reliable calculations up to a seawater ionic strength of 0.7 (Ball and Nordstrom 1991), it is inappropriate to apply geochemical modeling to salty groundwaters due to the poor estimation of the mixing proportion of seawater in cases with high EC (Fig. 5). In this study, the models proceeded only for groundwaters with seawater mixing proportions of <80%. Accordingly, the aquifer connectivity along four sections can be constructed, as shown in Fig.7, which provides complementary information on the devel-opment of a conceptual model of groundwater ﬂow.
Based on the results of geochemical modeling, the NCPGD can be divided into three aquifers. As shown in Fig.7, the spatial correlation of aquifers is roughly identical to chronological stratigraphy. The top aquifer is composed of Holocene sediments, the intermediate aquifer was deposited between 30,000 and 10,000 years ago (a), and the bottom of the deep aquifer lies on a temporal surface of 50,000 years ago (a). Subsequently, aquifer correlations (lines with question marks in Fig. 7) can be extended if geochemical modeling is not available. Aquifer thickness generally increases southward. In sections III and IV (Figs 7c and
7d, respectively), the bottom of the deep aquifer is too deep to appear, and the intermediate aquifer can be further divided into two subsystems in the downgradient area, whereas it remains unlayered in the upgradient area (Fig.7).
Chronostratigraphy and aquifer correlation
The aquifer is traditionally deﬁned as a composite of highly permeable layers. Therefore, the establishment of a conceptual groundwater model substantially relies on lithological correlation among wells. In this study, the results of geochemical modeling demonstrate that the hydraulic connectivity of an aquifer is highly dependent on chronological stratigraphy, which can be explained with sequence stratigraphy, as shown in Fig.8. The grain size of sediments in the coastal area is principally controlled by eustatic sea level and tectonic movement. The overall combination effect is that the depositional environment changes both in time and space. Two sand bodies deposited during different low-sea-level periods (marked with solid circles in Fig.8) are very likely to be separated by a clay layer during a high-sea-level period. Although the marked sand bodies are located at similar depth, they lie on different time horizons and are not interlinked in the groundwaterﬂow system. It is suggested that aquifers should be connected mostly based on the concept of chronological stratigraphy. The faunal/ﬂoral content and dating of sediments may serve as useful time markers to determine the aquifer correlation, especially for low-permeable, unconsolidated aquifers such as those of the NCPGD.
Estimation of aquifer hydraulic conductivity
For groundwaterﬂow on a regional scale, the radiocarbon technique is traditionally employed to reveal aquifer hydraulic conditions, with mean residence time of recharged water based on Darcy’s law (e.g., Cronin et al.
2005; Bouhlassa and Aiachi 2002; Buckau et al. 2000). Ignoring the effects of hydrodynamic dispersion, radio-carbon concentration in groundwater is mainly affected by radioactive decay and mass transfer between soil and groundwater (e.g., Castro and Goblet2005). It is obvious that a reliable geochemical reaction model between two monitoring wells is essential to calculate mean travel time. As listed inAppendix, some geochemical models involve large amounts of mass transfer, which may introduce error in the estimation of travel time. In addition, hydraulic heads in the NCPGD are highly variable both in time and space due to extensive groundwater withdraw. Under the circumstances, aquifer hydraulic conductivity can be only roughly estimated.
According to the results of inverse geochemical modeling as listed in Appendix, the ﬂux in groundwater can be derived from dividing distance between monitoring wells by modeled travel time. The annual mean of hydraulic gradient in the sampling year is also available from the Inquiry and Supply System of Water Resources Database (Water Resources Agency 2006). Applying Darcy’s law to aquifers along sections I–IV (Fig. 1), aquifer hydraulic conductivities can be subsequently derived (Table 4). As shown in Fig. 9, the intermediate aquifer in the NCPGD has been seriously overpumped, generating an extensive depression cone, especially around the areas of sites N07, N11, and N21 (Fig. 9). Table 3 Presumed chemistry of seawater and inﬁltration water
Items Inﬁltration water Seawater Temperature (°C) 25.0 25.0 pH 8.3 8.5 Na+(mg/l) 0.0 10,500.0 K+(mg/l) 0.0 390.0 Ca2+(mg/l) 19.7 410.0 Mg2+(mg/l) 0.0 1,350.0 Alkalinity (mg/l) 59.8 142.0 Cl−(mg/l) 0.0 19,000.0 SO42−(mg/l) 0.0 2,700.0 δ13C (‰) −12.5 0.0 δD (‰) −35.4 0.0 δ18 O (‰) −6.4 0.0 14 C (pmc) 75.0 14C of sediment
Fig. 7 Aquifer correlation along four sections, I (a) to IV (d), in the North Chianan Plain Groundwater District (NCPGD) (Fig. 1). The lines indicate the bottom limits for aquifers. The lines with question marks show the correlation derived from the concept of chronological stratigraphy
Thus, three segments around sites N11 and N21 have negative hydraulic gradients (Table4). It is also conﬁrmed that their upstream segments (N12-1→N11-2 and N12-2→N11-2) have hydraulic conductivities of nearly one order of magnitude higher than the others. In addition, the results demonstrate that the top aquifer is slightly more permeable than the intermediate and deep aquifers. The continuous recharge to the top aquifer is considered to be the source of error, whereas hydraulic conductivities are calculated under constant ﬂux. The extra recharge along the ﬂow path not only generates a more gentle hydraulic gradient but also supplies modern radiocarbon to reduce the apparent travel time. Both effects will contribute to higher hydraulic conductivities. On the other hand, the top aquifer can be deﬁned as an unconﬁned aquifer according to the higher derived hydraulic conductivities of the top aquifer.
To conclude the above discussion, the calculated values of 3.14× 10−5–1.80×10−4 m/s represent the range of hydraulic conductivities for unconsolidatedﬁne-sand aqui-fers in the NCPGD (Table 4). The US Army Corps of Engineers (1999) presented considerably lower values of 10−7–10−5 m/s for a common range ofﬁne-sand aquifers. Furthermore, an in situ single-well pumping test after well development and performance testing (Cooper and Jacob
1946; Taiwan Sugar Corporation1999) also obtained much smaller hydraulic conductivities (Table4). It is apparent that reverse geochemical modeling overestimates hydraulic con-ductivities by roughly one order of magnitude. Generally, an in situ pumping test provides genuine properties of aquifer permeability, whereas radiocarbon aging gives a combined result of all inﬂuences on groundwater ﬂow, such as aquifer permeability and artiﬁcial extraction of groundwater. In the Fig. 8 Sketch of conceptual sequence stratigraphy (Van Wagoner et al. 1990). The black dots represent sand bodies deposited during different low-sea-level periods
Table 4 Calculation of aquifer hydraulic conductivity
Well pairs Aquifer Da(km) Travel time (years) Ib(×10−3) K (m/s)ccalculated K (m/s)din situ testing N03-1→N02-2 Intermediate 10.7 5,447 1.98 3.14×10−5 6.97×10−5 N03-2→N02-2 Intermediate 11.0 1,803 1.68 1.15×10−4 6.97×10−5 N02-1→N01-1 Top 10.7 764 0.99 4.48×10−4 – N09-1→N08-2 Intermediate 10.1 2,949 2.81 3.86×10−5 2.08×10−5 N09-3→N08-3 Deep 9.7 3,195 1.47 6.54×10−5 2.16×10−5 N08-1→N07-1 Top 6.7 5,369 0.23 1.72×10−4 6.05×10−5 N08-2→N07-3 Intermediate 7.5 1,497 0.88 1.80×10−4 2.13×10−5 N12-1→N11-2 Intermediate 8.5 428 0.52 1.21×10−3 1.95×10−6 N12-2→N11-2 Intermediate 8.0 424 0.56 1.07×10−3 1.95×10−6 N12-2→N11-3 Intermediate 9.3 899 2.46 1.33×10−4 3.12×10−5 N11-1→N10-1 Top 13.2 5,076 0.10 8.20×10−4 4.84×10−5 N11-2→N10-3 Intermediate 13.3 3,559 Negative – 4.96×10−5 N11-3→N10-4 Intermediate 13.6 5,936 Negative – 5.57×10−5 N10-4→N13-1 Intermediate 10.4 5,733 0.37 1.55×10−4 1.88×10−5 N10-4→N13-2 Intermediate 10.5 5,495 0.52 1.16×10−4 2.22×10−6 N22-1→N21-2 Intermediate 5.9 351 5.06 1.05×10−4 4.70×10−5 N21-2→N20-2 Intermediate 12.5 9,402 Negative – 1.63×10−4 a
D (km), represents distance between well pair b
I (×10−3), represents hydraulic gradient between well pair c
K (m/s) calculated, represents the hydraulic conductivity derived by reverse geochemical modeling d
NCPGD, the substantial extraction of groundwater over a long period could accelerate groundwaterﬂow, leading to an overestimation of aquifer permeability.
The NCPGD is principally composed of low-permeable aquifers, and the natural inﬁltration of fresh groundwater is unable toﬂush the aquifers. The results of geochemical analysis demonstrate that the aquifers contain considerable amounts of old brine water. The high content of sulfur causes a highly reduced environment to release arsenic from muddy sediments. However, the high variation of the geochemical environment makes geochemical modeling practicable to track the groundwaterﬂow path.
Although the stable isotopic compositions of oxygen and hydrogen are not available, the seawater mixing proportion can be estimated using the relationship between mixing proportion of seawater and EC derived from the isotopic compositions of the SCPGD. The resultant mixing
propor-tion was regarded as a constraint in the NETPATH model. The hydraulic connectivity was examined between pairs of possible interlinked wells for four sections along the groundwater ﬂow path, and the aquifer conceptual model was successfully constructed. It is suggested that aquifer correlation is mostly based on chronological stratigraphy, especially for the low-permeable, unconsolidated aquifers.
Darcy’s law was applied along the direction of ground-water ﬂow, and the hydraulic conductivities of NCPGD aquifers were roughly estimated to be in the range of 3.14× 10−5–1.80×10−4 m/s. Compared with in situ single-well pumping tests after well development and performance test, the reverse geochemical modeling tends to overestimate aquifer hydraulic conductivity by roughly one order of magnitude. It could be the result of extensive groundwater withdrawal over a long period in the NCPGD.
Acknowledgements Special thanks go to Mr. L.Y. Fei, who promoted this project. The authors appreciate helpful comments from reviewers and journal editors. This study was supported by research grants from Water Resources Agency and Central Geological Survey, Ministry of Economic Affairs, Taiwan. Fig. 9 Contour map showing annual mean of hydraulic head (in meter above sea level) in the sampling year. The contour interval is 5 m. Creeks and lakes are indicated (see Fig.1for labels)
Results of NETPATH models
Model (pairs of wells) Isotopes Mixing proportion
δ13 C (‰) Observed δ13 C (‰) Travel Time (years)
Source Seawater Inﬁltration
N03-1→N02-1 No −5.8 −10.0 −757 0.289 0.000 0.711 N03-1→N02-2 Yes −17.8 −18.2 5,447 0.990 0.010 0.000 N03-2→N02-2 Yes −18.6 −18.2 1,803 0.990 0.010 0.000 N02-1→N01-1 Yes -11.1 −10.9 764 0.996 0.004 0.000 N02-2→N01-2 No −8.8 −12.9 4,498 0.989 0.011 0.000 N09-1→N08-1 No No model obtained N09-1→N08-2 Yes 5.9 7.3 2,949 0.953 0.047 0.000 N09-2→N08-2 No 1.2 7.3 3,083 0.980 0.020 0.000 N09-2→N08-3 No −15.5 −13.7 −63 0.948 0.052 0.000 N09-3→N08-4 No −3.9 −6.2 −16,954 0.908 0.092 0.000 N08-1→N07-1 Yes −23.5 −24.8 5,369 0.170 0.830 0.000 N08-2→N07-2 Yes 6.5 6.2 1,257 0.867 0.133 0.000 N08-2→N07-3 Yes 6.0 5.4 1,497 0.973 0.027 0.000 N08-3→N07-3 No 0.5 5.4 −11,139 0.942 0.058 0.000 N08-4→N07-4 No −22.4 −16.5 −749 0.993 0.007 0.000 N07-3→N06-1 No No model obtained N07-4→N06-1 Yes −11.9 −11.0 5,099 0.987 0.013 0.000 N12-1→N11-1 No −8.9 −8.8 −4,015 0.696 0.304 0.000 N12-1→N11-2 Yes 5.2 4.8 428 0.718 0.282 0.000 N12-2→N11-2 Yes 5.5 4.8 424 0.720 0.280 0.000 N12-2→N11-3 Yes −13.9 −16.0 899 0.955 0.045 0.000 N12-3→N11-3 No −17.8 −16.0 −10,129 0.939 0.061 0.000 N11-1→N10-1 Yes −22.5 −21.6 5,076 0.122 0.878 0.000 N11-2→N10-2 No −1.5 −0.2 −385 0.452 0.548 0.000 N11-2→N10-3 Yes 6.6 9.9 3,559 0.992 0.008 0.000 N11-3→N10-4 Yes −19.3 −17.5 5,936 0.995 0.005 0.000 N10-3→N13-1 No No model obtained N10-4→N13-1 Yes −6.2 −5.9 5,733 0.866 0.134 0.000 N10-4→N13-2 Yes −10.0 −10.3 5,495 1.000 0.000 0.000 N22-1→N21-1 No −19.8 −8.6 −2,172 0.984 0.016 0.000 N22-1→N21-2 Yes −3.0 −2.3 351 0.995 0.005 0.000 N21-2→N20-2 Yes −1.5 −0.8 9,402 0.858 0.142 0.000 N21-3→N20-2 No No model obtained N21-3→N20-3 No −7.6 −6.9 −1,174 0.864 0.136 0.000 N20-2→N19-3 No No model obtained N20-3→N19-3 No No model obtained
Mass transfer of phases (mmol/l)
(no sign): phase precipitation from solution; (-): phase dissolution into solution
Calcite Dolomite Goethite Pyrite Exchange Gypsum Siderite CH2O CH4 CO2 H2S
0.000 1.350 0.000 0.000 0.696 0.000 0.000 0.000 0.188 0.000 0.000 0.000 −0.775 −0.070 0.000 3.467 2.630 0.000 0.000 2.874 −0.107 −2.885 0.000 −0.131 0.018 0.000 3.526 2.992 0.000 6.498 0.000 −0.975 −3.242 1.174 0.099 0.057 0.006 0.801 1.352 0.000 0.000 0.000 1.574 0.000 3.313 0.692 0.067 0.000 2.785 0.000 0.000 0.000 0.370 −6.399 −0.360 No model obtained 4.486 −1.568 0.000 0.000 3.272 0.000 0.000 4.456 −3.030 0.000 0.000 12.788 0.000 0.000 0.000 13.051 0.000 0.000 0.000 0.000 -8.320 0.000 1.680 −2.403 0.000 0.000 0.000 0.000 0.000 0.000 −0.006 6.319 0.000 0.000 −4.052 0.000 0.000 −3.078 0.000 0.000 26.490 −13.253 0.000 0.000 27.439 −19.149 6.900 0.000 10.806 0.000 −6.885 0.000 22.379 0.000 −21.515 0.992 −1.463 0.000 0.071 −6.659 −3.989 0.000 0.000 −0.121 0.000 0.000 0.000 0.325 4.115 0.000 −3.714 −1.282 −4.059 0.000 0.000 0.000 0.520 −1.638 −0.793 −11.677 −0.831 −5.054 0.000 12.566 0.000 0.000 −11.068 0.000 −15.431 −2.362 0.188 0.000 −16.948 0.000 0.000 0.000 0.032 19.262 0.000 No model obtained 1.029 −0.396 −0.201 0.000 0.000 −0.500 0.000 0.000 −0.027 −3.962 0.000 −1.163 0.963 0.275 0.000 −1.286 0.567 0.000 0.079 0.000 0.000 0.000 47.704 −4.999 0.000 0.036 34.327 −8.230 0.000 0.000 0.000 −36.629 −0.055 46.800 −4.591 0.041 0.001 33.313 −8.209 0.000 0.000 0.000 −37.771 0.000 1.718 −0.059 0.000 0.038 0.000 0.000 0.000 0.000 1.379 −4.567 −1.448 −1.144 0.036 0.000 −3.062 −0.773 0.000 2.456 0.000 0.000 −1.228 0.000 21.387 −16.967 0.000 0.027 4.173 0.000 0.000 0.000 15.830 0.000 −15.879 20.294 −10.498 0.000 −0.004 −7.987 −15.758 0.063 0.000 0.000 0.000 0.000 0.000 5.830 0.000 0.000 3.877 0.000 0.017 0.000 0.000 −11.939 −0.007 0.587 −2.841 −1.015 0.000 0.000 0.000 0.996 0.000 0.000 7.463 −0.140 No model obtained 4.476 −1.680 0.000 0.000 2.039 −1.229 −0.014 16.264 −8.048 0.000 0.000 7.607 −0.122 0.241 0.021 7.641 0.000 0.000 0.000 0.000 −0.743 0.000 −0.417 −0.984 0.000 0.000 0.000 0.000 0.000 0.000 0.000 7.384 0.000 6.732 0.000 0.007 0.000 7.596 −0.149 0.000 6.713 −3.348 0.000 0.000 2.803 −0.077 −0.014 0.000 0.000 −2.267 0.000 0.000 0.000 −0.888 0.154 No model obtained 5.556 −0.062 0.000 0.011 3.351 −2.167 0.000 0.000 0.000 0.000 0.130 No model obtained No model obtained
André L, Franceschi M, Pouchan P, Atteia O (2005) Using geochemical data and modelling to enhance the understanding of groundwaterﬂow in a regional deep aquifer, Aquitaine Basin, south-west of France. J Hydrol 305:40–62
Appelo CAJ (1994) Cation and proton exchange, pH variations, and carbonate reactions in a freshening aquifer. Water Resour Res 30(10):2793–2805
Bajjali W (2006) Recharge mechanism and hydrochemistry evalu-ation of groundwater in the Nuaimeh area, Jordan, using environmental isotope techniques. Hydrogeol J 14:180–191 Ball JW, Nordstrom DK (1991) User’s manual forWATEQ4F, with
revised thermodynamic data base and test cases for calculating speciation of major, trace, and redox elements in natural waters. U.S. Geological Survey Open-File Report 91–183
Bouhlassa S, Aiachi A (2002) Groundwater dating with radiocarbon: application to an aquifer under semi-arid conditions in the South of Morocco (Guelmime). Appl Radiat Isotopes 56:637–647 Buckau G, Artinger R, Geyer S, Wolf M, Fritz P, Kim JI (2000)14C
dating of Gorleben groundwater. Appl Geochem 15:583–597 Castro MC, Goblet P (2005) Calculation of ground water ages– A
comparative analysis. Ground Water 43(3):368–380
Central Geological Survey (2006) Hydrogeology Database. http:// hydro.moeacgs.gov.tw. Cited 4 Aug 2006
Chen WS, Yang CC, Yang HC, Wu LC, Lin CW, Chang HC, Shih RC, Lin WH, Lee YH, Shih TS, Lu SD (2004) Tectono-geomorphologic studies in the Chiayi-Tainan region in south-western Taiwan and its implicaitons for active structure. Bulletin of the Central Geological Survey, Taiwan 17:53–77 [in Chinese]
Coleman ML, Shepherd TJ, Durham JJ, Rouse JE, Moore GR (1982) Reduction of water with zinc for hydrogen isotope analysis. Anal Chem 54(6):993–995
Cooper HH, Jacob CE (1946) A generalized graphical method for evaluating formation constants and summarizing well ﬁeld history. Am. Geophys. Union Trans. 27:526–534
Cronin AA, Barth JAC, Elliot T, Kalin RM (2005) Recharge velocity and geochemical evolution for the Permo-Triassic Sherwood Sandstone, Northern Ireland. J Hydrol 315:308–324 Epstein S, Mayeda T (1953) Variation of O-18 content of waters
from natural sources. Geochim Cosmochim Acta 4:213–224 Goldberg ED, Broecker WS, Gross MG, Turekian KK (1971)
Marine chemistry. In: Radioactivity in the marine environment. Washington, D.C., National Academy of Sciences, pp 137–146 Güler C, Thyne GD (2004) Hydrologic and geologic factors controlling surface and groundwater chemistry in Indian Wells-Owens Valley area, southeastern California, USA. J Hydrol 285:177–198
Hiatt EE, Kyser K, Dalrymple RW (2003) Relationships among sedimentology, stratigraphy, and diagenesis in the Proterozoic Thelon Basin, Nunavut, Canada: implications for paleoaquifers and sedimentary-hosted mineral deposits. J Geochem Explor 80:221–240
Hidalgo MC, Cruz-Sanjulián J (2001) Groundwater composition, hydrochemical evolution and mass transfer in a regional detrital aquifer (Baza basin, southern Spain). Appl Geochem 16:745–758
Ingerson E, Pearson Jr. FJ (1964) Estimation of age and rate of motion of groundwater by the 14C-method. In: Recent researches in theﬁelds of atmosphere, hydrosphere, and nuclear geochemistry, Sugawara Festival Volume. Maruzen Co., Tokyo, pp 263–283
Lu FJ (1990) Blackfoot disease: arsenic or humic acid ? Lancet 336:115–116
Lubczynski M, Roy J (2003) Hydrogeological interpretation and potential of the new magnetic resonance sounding (MRS) method. J Hydrol 283:19–40
Parkhurst DL, Appelo CAJ (1999) User’s guide to PHREEQC (version 2)– A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Water-Resources Investigation Report 99-4259, Denver
Peng TR (2001) Study of stable hydrogen and oxygen isotopes for groundwaters in Chianan Plain and Ilan Plain, Taiwan Ground-water Monitoring Network Project. Central Geological Survey, Ministry of Economic Affairs, Taiwan [in Chinese]
Plummer LN, Prestemon EC, Parkhurst DL (1994) An interactive code (NETPATH) for modeling net geochemical reactions along a ﬂow path, version 2.0. U.S. Geological Survey Water Resources Investigations Report 94-4169, Reston
Shtivelman V, Goldman M (2000) Integration of shallow reﬂection seismics and time domain electromagnetics for detailed study of the coastal aquifer in the Nitzanim area of Israel. J Appl Geophys 44:197–215
Sugarman PJ, Miller KG (1997) Correlation of Miocene sequences and hydogeologic units, New Jersey Coastal Plain. Sediment Geol 108:3–18
Taiwan Sugar Corporation (1999) Construction of groundwater monitoring well and related hydraulic tests, Taiwan Groundwa-ter Monitoring Network Project. WaGroundwa-ter Resources Agency, Ministry of Economic Affairs, Taiwan [in Chinese]
US Army Corps of Engineers (1999) Occurrence and movement of groundwater. In: Groundwater hydrology. Engineer Manual 1110-2-1421
Van der Kemp WJM, Appelo CAJ, Walraevens K (2000) Inverse chemical modeling and radiocarbon dating of palaeoground-waters- the Tertiary Ledo-Paniselian aquifer in Flanders, Belgium. Water Resour Res 35(3):1277–1287
Van Geen A, Zheng Y, Versteeg R, Stute M, Horneman A,Dhar R, Steckler M, Gelman A, Small C, Ahsan H, Graziano JH, Hussain I, Ahmed KM (2003) Spatial variability of arsenic in 6000 tube wells in a 25 km2area of Bangladesh. Water Resour Res 39(5):1140–1155
Van Wagoner JC, Mitchum RM, Campion KM, Rahmanian VD (1990) Siliciclastic sequence stratigraphy in well logs, cores, and outcrops. American Association of Petroleum Geologists Methods in Exploration Series 7, Tulsa
Water Resources Agency (2006) Inquiry and Supply System of Water Resources Database.http://gweb.wra.gov.tw. Cited 4 Aug 2006 Wigley TML, Plummer LN, Pearson FJ (1978) Mass transfer and
carbon isotope evolution in natural water systems. Geochim Cosmochim Acta 42:1117–1139
Zuber A, Weise SM, OsenbruÈck K, Pajnowska H, Grabczak J (2000) Age and recharge pattern of water in the Oligocene of the Mazovian basin (Poland) as indicated by environmental tracers. J Hydrol 233:174–188