Spatial Variability of Groundwater Arsenic Concentration as Controlled by Hydrogeology: Conceptual Analysis Using 2-D Reactive Transport Modeling.

Combined geological, hydrogeological, and geochemical controls on the arsenic concentration of contaminated aquifers in SE Asia were explored by two-dimensional (2-D) reactive transport modeling of data sets from Bangladesh, Cambodia, and Vietnam. For each site, the field data are summarized and used to create a conceptual 2-D reactive transport model that elucidates characteristic features influencing the groundwater arsenic concentration. Comparison of models for Bangladesh and Vietnam indicates that fine-grained layers overlying young sandy aquifers generate shallow high arsenic groundwater because low vertical groundwater velocities allow sufficient time for kinetic As release from the sediment. The low vertical groundwater velocity below major river channels, predicted by the model, also creates long groundwater residence times, leading to high arsenic groundwater. Young aquifer sediments release more arsenic than older sediments, and alternating young and older sediments create complex patterns of high and low arsenic groundwater. Over time, floodplain basins will subside, and river channels migrate, causing sedimentation and erosion on the floodplain while creating local environments with evolving hydrogeology and groundwater geochemistry. We have developed a three-step model for the evolution of the Red River floodplain with sedimentation and shifting channels over the last 6000 years. The results show comparable timescales between the dynamics of arsenic release and of river migration, causing complex groundwater As distributions, comprising geochemical palinopsia of long vanished rivers.

Introduction

The high spatial variability of the groundwater arsenic content over even small distances has been an issue of concern because it makes the management of the groundwater resource very difficult (British Geological Survey/Department of Public Health Engineering, 2001; van Geen et al., 2003; Winkel et al., 2011). The variability of the groundwater arsenic concentration depends, among other factors, on the sedimentary environment (Acharyya et al., 2000; Donselaar et al., 2017; Hoque et al., 2012; Papacostas et al., 2008; Sahu & Saha, 2015; Weinman et al., 2008) and the burial age of the sediments, which affect the available amounts, and reactivities, of organic matter and Fe oxides (Kocar et al., 2008; Postma et al., 2012; Radloff et al., 2017; Stuckey et al., 2015). In addition, hydrogeological conditions, such as recharge and hydraulic conductivity, affect the spatial distribution of recharge rates and flow patterns, as well as the number of pore water volumes flushed through the deposits. This flushing drives the leaching of adsorbed arsenic (Desbarats et al., 2014; Postma et al., 2016; Radloff et al., 2017; Sø, Postma, Vi, Pham, Kazmierczak, et al., 2018; van Geen et al., 2008; Weinman et al., 2008). The spatial variation of all these interacting variables gives rise to the observed highly complex As distribution patterns and large variability of concentrations even on a local scale. Reactive transport modeling lends itself to the evaluation of integrated effects of variations in geology, hydrogeology, and geochemical processes on the groundwater arsenic distribution. Comprehensive reactive transport modeling that includes all the main geochemical processes has only been applied in a few studies of As mobility (Kocar et al., 2014; Postma et al., 2007, 2016; Rawson et al., 2017). A prerequisite for performing meaningful comprehensive reactive transport modeling is the availability of data sets that comprise both the detailed geology and burial age of the deposits, groundwater flow patterns, and groundwater ages as well as a detailed description of the groundwater chemistry. We have searched the literature for such data sets and found several well‐documented hydrogeological settings that yield a distinct groundwater arsenic distribution. We use two‐dimensional (2‐D) reactive transport modeling to analyze the effects of hydrogeological parameters on the spatial distribution of the arsenic concentration and try to identify characteristic patterns belonging to a given hydrogeological setting. We also demonstrate how the groundwater arsenic distributions at the different field sites can be quantified by using the same 2‐D reactive transport model that is based on the one‐dimensional (1‐D) model of Postma et al. (2016). These cases provide conceptual elements for the hydrogeological control on groundwater arsenic distributions. The elements can be recognized elsewhere on floodplains and used as building blocks for comprehending patterns of groundwater arsenic distribution at regional scales. Once the model has been validated on current hydrogeological settings, we use it for modeling the groundwater As distribution during the geological history of an evolving floodplain. Toward this purpose, we have constructed a tentative model for the evolution of the groundwater arsenic distribution on the Red River floodplain as a function of geological and hydrogeological changes over the last 6000 years.

Model Development

The model used in this study is PHAST‐3, which combines the flow model HST3D (without its density modeling features) and the general geochemical model PHREEQC. PHAST‐3 has the same functionality as PHAST‐2 (Parkhurst et al., 2010) but makes use of the latest PHREEQC‐3 (Parkhurst & Appelo, 2013) based PhreeqcRM reaction module. The geochemical model used in PHAST‐3 is the same as the 1‐D PHREEQC model described in Postma et al. (2016) for aquifers on the Red River floodplain, Vietnam. The only difference in the PHAST model, compared to the 1‐D PHREEQC model, is that the more complex geology in the final 2‐D model includes Pleistocene clay for which the initial content of organic carbon and Fe oxides was divided by 10 and Pleistocene sand where the initial values were divided by 100. While the Pleistocene sediments generally are recognized to be oxidized and therefore less reactive, the division numbers used here are arbitrary. Like the 1‐D model, dispersion has not been included. Completely removing numerical dispersion from the 2‐D models would result in unmanageable long computational runs, and it was decided not to add additional dispersion.

The geochemical model used contains a description of how the reaction rates of organic matter degradation and reduction of arsenic‐bearing Fe oxides change over time. The model was calibrated on the water chemistry in depth profiles of three aquifers with a burial age of, respectively, 515, 3530, and 5900 years on the Red River floodplain, Vietnam (Postma et al., 2016). Here it is applied unaltered, using the same rate parameters and sediment concentrations, to model field data from the various field sites in Bangladesh, Cambodia, and Vietnam. Below is a brief summary of the geochemical model highlighting the processes required to understand the model results presented here. For a more complete description, see Postma et al. (2016).

The main process mobilizing arsenic from the sediment into the groundwater is the reductive dissolution of arsenic containing Fe oxide, mainly goethite:

Reaction (1) couples the oxidation of organic carbon to the reduction of Fe oxides while releasing Fe(II) and trace amounts of arsenic contained in the Fe oxides. The rate of organic carbon degradation, the overall reaction control, is given by the following equation:

Here m  = 1.36 mol/L is the initial amount of organic matter in the sediment, based on the analytical concentration recalculated to the concentration per liter of contacting groundwater, m is the remaining amount of organic carbon in the sediment at time t, and 9.3 × 10−12 s−1 is the rate constant. The exponent to (m /m ) reflects the variation width of the organic carbon reactivity. For a first‐order rate dependency, the exponent to (m /m ) would be 1. The exponent value of 2.5 reflects a larger difference between the most and least reactive organic carbon; that is, the initial pool of organic carbon is more reactive than the remaining part at a later time. The amount of Fe oxides that reduces is given by

Here m , the initial mass of FeOOH in the sediment, is 0.394 mol/L corresponding to 65 μmol FeOOH/g sediment, obtained by extrapolation from the field data, and m is the remaining amount of FeOOH in the sediment at time t. SR is the saturation ratio, the ion activity product (for the reaction FeOOH + 3H+ ↔ Fe3+ + 2H2O) divided by the equilibrium constant for the reaction. SR depends on the redox conditions, controlling the Fe3+ activity, and the pH with a low pH favoring FeOOH reduction. The pH is again a function of the reduction of FeOOH itself, the precipitation of FeCO3, and the dissolution of CaCO3. The rate of organic carbon oxidation is always larger than the rate of FeOOH reduction, and the excess electron flow from organic carbon will, in these low sulfate aquifers, appear in the groundwater as dissolved CH4 following the overall reaction:

The Fe oxide that is reduced contains 1.2 mmol As(V) per mol of Fe and is the source of As(III) released to the groundwater. Once in the groundwater, As(III) may adsorb to the sediment as described by the Langmuir adsorption model of Nguyen et al. (2014). The Langmuir equation states that

Here reflects the sorbed As(III), [S tot] the total number of available sites, K = 103.176 the adsorption constant as determined by Nguyen et al. (2014), and [H3AsO4] the activity of H3AsO4. Experiments have shown [S tot] to be related to the Fe oxide content of the sediment (Sø, Postma, Vi, Pham, Viet, & Jakobsen, 2018). The number of adsorption sites is therefore coupled to the amount of Fe oxide present with a site density of 0.07 moles of adsorption sites per mole of Fe oxide. Accordingly, the number of adsorption sites decreases as more Fe oxide is reduced. Reactive transport modeling employed the WATEQ4F.DAT database with the compilation of Langmuir et al. (2006) for arsenic speciation. The ordinary differential equation solver (CVODE) for stiff systems was used for the kinetic calculations.

The coupled kinetic reactions slow down the geochemical calculations, but PHAST can run in Message Passing Interface parallel programming. PHAST uses operator splitting, solving the transport and the chemistry equations sequentially without iteration between the two, making parallelization efficient. Parallelization is optimized by reallocation of processing units to where they are most needed after each time step. When the number of parallel processes exceeds 20, the reduction in calculation time per added parallel process decreases, placing limits on how large a model can be solved within a reasonable time. Model steps used in this study required up to 60 hr using 40 hyperthreaded cores running at 2.3 GHz.

PHAST‐3 can transfer a set of geochemical parameters from one model to another by starting a model with the geochemical results from a preceding model as the initial state. The geometry and the grid of the new model can be different, because geochemical parameters are only transferred, by interpolation, to the part of the model volume that is retained in the subsequent model. The evolution of a sedimentary sequence can thus be modeled stepwise, where the sediment and water chemistry from the volume that is not eroded are inherited as the initial conditions for the following step in the hydrogeological evolution.

Sø, Postma, Vi, Pham, Kazmierczak, et al. (2018) introduced the concept of total numbers of pore volumes flushed to describe the leaching of As from the sediment. The idea behind the concept is that As‐free water is recharged from the surface and flushes As from the sediment. The analysis of Sø, Postma, Vi, Pham, Kazmierczak, et al. (2018) only considers the vertical flow component. However, in a 2‐D model it does not correctly describe the total number of times the groundwater has actually been replaced in the porespace. Therefore, we here rephrase the concept as the total number of recharged pore volumes (RPVs). From field data, the total number of RPVs flushed through the aquifer since deposition is calculated by multiplying sediment burial age (year) with groundwater velocity (m/year) and dividing by the traveled vertical distance (m) (Sø, Postma, Vi, Pham, Kazmierczak, et al., 2018, equation 9), and it is equal to the ratio of sediment burial age divided by the groundwater age at a given depth. To keep track of the RPVs, a kinetic function was written within PHREEQC. For every year, the function adds 1E−8 to the concentration of a defined component named Age and another component named Sda (sediment age). Age accumulates in the groundwater over time giving the groundwater age, while Sda precipitates as a mineral phase named SedAge with a very low solubility. The number of RPVs is found by dividing SedAge with Age.

Results

Here we present a summary of the geochemical and hydrological properties of four well‐studied field sites in SE Asia, followed by a conceptual 2‐D reactive transport model that emphasizes the features we find characteristic and of general interest. The objective is to elucidate the hydrogeological parameters that influence the spatial distribution of arsenic in the groundwater. For reasons of space and clarity, we have confined ourselves to displaying only the arsenic model results even though the model produces results for all bulk chemical parameters.

H‐Transect at DanPhuong, Red River Floodplain, Vietnam

Field Data

The field site is located 30 km upstream from Hanoi (105°38′06″E, 21°09′01″ N) on the Red River floodplain between the river and the dyke where the land surface is cultivated without paddy fields (Larsen et al., 2008; Postma et al., 2007). Mean annual precipitation is 1,640 mm with 70% falling in the monsoon between May and October and with a potential annual evaporation of 700–800 mm (Larsen et al., 2008). The sandy aquifer has a burial age between 460 and 600 years, as determined by optical stimulated luminescence (Sø, Postma, Vi, Pham, Kazmierczak, et al., 2018), and it is overlain by a highly fractured clay layer. Dating by 3H/3He of the groundwater in the aquifer indicates a downward vertical groundwater velocity of 0.5 m/year (Postma et al., 2007, 2012). Throughout this paper, groundwater velocity is the actual velocity of water through the pores. The total number of RPVs is between 22 and 39 (Sø, Postma, Vi, Pham, Kazmierczak, et al., 2018). The idea behind using RPVs is that they reflect the number of pore volumes, initially recharged as arsenic‐free water, that pass through the deposit. The water chemistry displays the effects of the degradation of organic carbon associated with the reductive dissolution of As‐bearing Fe oxides as increases in As, Fe(II), and CH4 with depth (Postma et al., 2007). The As concentration (Figure 1) reaches a maximum of 6 μmol/L at about 10 m below the water table, and arsenic, as well as the other parameters, show a distinct horizontal layering in the concentration distribution with generally increasing concentrations over depth (Postma et al., 2007).

Figure 1

Groundwater arsenic in a cross section through the H‐transect, a sandy aquifer on the shore of the Red River at Dan Phuong, Vietnam. Crosses indicate sampling points (based on Postma et al., 2007).

Modeling

Postma et al. (2016) developed a 1‐D reactive transport model for aquifers of different burial ages on the Red River floodplain, which was also applied to the H‐transect. It is a far more comprehensive model than the 1‐D reactive transport model developed only for the H‐transect by Postma et al. (2007). The geochemical model, developed in the 1‐D model of Postma et al. (2016), is used unmodified in the 2‐D models presented here. A 2‐D PHAST model, covering the H‐transect, was set up with 3,600 m in total length, the distance from the water divide to the river, and a depth of 45 m derived from Larsen et al. (2008). The model section displayed in Figure 2a is the distance interval 3,400–3,500 m, to a depth of 20 m, corresponding to the field section presented in Figure 1. A homogenous recharge rate of 175 mm/year is imposed over the whole aquifer surface and with a porosity of 0.35; this corresponds to a vertical downward groundwater velocity of 0.5 m/year in agreement with 3H/3He dating of the groundwater (Postma et al., 2007, 2012). The model was run for 530 years corresponding to the time from a deposition of the entire deposit at the average burial age of the aquifer sediment.

Figure 2

Conceptual two‐dimensional (2‐D) reactive transport model for the H‐transect, Vietnam (Figure 1). (a) Model setup, (b) As distribution, and (c) groundwater age and flow pattern. Groundwater flow, in this and subsequent plots, is depicted by vectors starting in a point, the line indicating flow direction and relative velocity by length. Shown is the distance interval 3,400–3,500 m of a 3,600 × 45‐m model with upstream no‐flow and downstream constant‐head boundary conditions.

The model results show a gradual increase in the groundwater As concentration over depth (Figure 2b), reaching about 4.5 μmol/L As(III) at 12‐m depth in groundwater 25 to 40 years old. This pattern is similar to that displayed by the field data (Figure 1). A striking common feature for both the field data and model results is the horizontal layering in the As concentration, which is a result of the age distribution of the groundwater (Figure 2c) that is again a function of the flow pattern. For a homogeneous sandy aquifer with uniform recharge, the groundwater will have the same age at each given depth with the groundwater age increasing over depth. This behavior can also be derived from a simple groundwater mass balance (Vogel, 1967; Appelo & Postma, 2005, p. 68). A direct result of the kinetic rate control on the reductive dissolution of As‐containing Fe oxides is that the model predicts the arsenic concentration to increase with increasing groundwater age. The high arsenic concentration building up in the H‐transect aquifer is the result of the young burial age of the aquifer sand, containing reactive sedimentary organic carbon and Fe oxide, and the limited rate of recharge through the fine‐grained top layer, implying a slow vertical groundwater velocity, allowing sufficient time for the kinetically released arsenic to build up. The low recharge rate at the same time delays the flushing of adsorbed As.

Araihazar Aquifer, Ganges‐Brahmaputra Floodplain, Bangladesh

The Araihazar aquifer is situated on the Ganges‐Brahmaputra floodplain, Bangladesh (90°37′43″ E, 23°47′50″N), and has been studied in great detail (Aziz et al., 2008; Horneman et al., 2004; Radloff et al., 2017; Stute et al., 2007; van Geen et al., 2003; Weinman et al., 2008]. On this floodplain, young sandy aquifer sediments, with a burial age of about 600 years, as determined by optical stimulated luminescence (Weinman et al., 2008), are partially covered by a fine‐grained toplayer. The groundwater arsenic concentration in this area ranges from <0.13 up to about 12 μmol/L As (Figure 3). Typical for this site is the high spatial variability in the groundwater arsenic concentration even over small distances. Weinman et al. (2008) and Aziz et al. (2008) found groundwater arsenic concentrations to be particularly low downstream of the sandy aquifer outcrops (Figure 3). They suggested that high recharge rates through the sandy windows could have caused flushing of the arsenic from the aquifer. In support, Stute et al. (2007) used 3H/3He dating of the groundwater to determine recharge rates of up to 1.1 m/year below sandy outcrops while recharge rates below fine‐grained layers are down to 0.1 m/year. They described a regional groundwater flow pattern toward the river based on data from 2003 to 2004. However, irrigation pumping may since have changed this general flow pattern (Radloff et al., 2017).

Figure 3

Groundwater arsenic in a cross section through the Araihazar aquifer, situated on the Ganges‐Brahmaputra floodplain, Bangladesh (based on Weinman et al., 2008).

A model was set up to schematically reflect the situation described by Weinman et al. (2008). A sandy aquifer with a burial age of 600 years is covered by a fine‐grained layer that contains several sand windows (Figure 4a). The recharge rates through the sand windows and the fine‐grained layer are, respectively, 1.1 and 0.1 m/year in accordance with the estimates of Stute et al. (2007). The model was run for 600 years, corresponding to the time from a deposition of the entire deposit at one time, to the burial age of the aquifer. The resulting arsenic distribution is depicted in Figure 4b. Model simulations show a low As concentration, <0.2 μmol/L, downstream of the sandy outcrop. This extends both directly below the sandy outcrop and below the upstream part of the clay/silt cover. Reversely, downstream of the clay/silt cover, high As concentrations can be found in the shallow sandy outcrop. The data resolution of the field data does not allow the recognition of the effects of these detailed flow patterns. Maximum modeled As concentrations are up to 1.6 μmol/L, which is lower than measured in the field that exceeds 10 μmol/L. The rate constants and amounts of organic matter and Fe oxide in the model are based on the Vietnam H‐transect, and the As release rate in the Araihazar aquifer is apparently higher. The modeled difference in groundwater As content downstream the two areas with a different recharge rate is the combined effect of groundwater transport and the geochemical processes. Downstream the fine‐grained layer, the low vertical groundwater flow rate of 0.4 m/year (porosity is set to 0.25) allows high As to build up in the groundwater (Figure 4b). In contrast, downstream the sand windows, the high vertical flow rates of 4.4 m/year dilute the As released from the sediment. The high vertical groundwater velocity also causes faster flushing of adsorbed arsenic from the sediment. Figure 4c shows the number of RPVs that have passed through the aquifer since burial 600 years ago. Below the sand windows, more than 400 RPVs have passed through the deposit, while below the fine‐grained layers, less than 100 RPVs have passed. The effect of different pore volumes is consistent with the findings of Sø, Postma, Vi, Pham, Kazmierczak, et al. (2018) who concluded that at least 200 RPVs are needed to flush arsenic from the Red River Delta aquifer.

Figure 4

Conceptual 2‐D reactive transport model of the Araihazar aquifer, Bangladesh. (a) Model setup, (b) As distribution, and (c) flow pattern and RPVs flushed since burial.

Tay Island, Mekong Floodplain, Vietnam

When reduced groundwater, enriched in As(III) and Fe(II), discharges into the oxidizing riverine environment, reactive As‐containing Fe oxides precipitate near the sediment surfaces (Berube et al., 2018; Datta et al., 2009). While groundwater discharge mainly occurs at the sides of the river, currents will redistribute the high As‐reactive sediment throughout the river bed, resulting in the occurrence of high As‐reactive river sediment in major rivers (Jung et al., 2012; Postma et al., 2010; Stahl et al., 2016). Within these large rivers, stationary or transient islands are present that will accumulate the high As‐reactive river sediments. Upon reduction, these sediments release As to the aqueous phase (Postma et al., 2010), as illustrated by reports of islands in major Asian rivers that contain high arsenic groundwater (Goswami et al., 2014; Papacostas et al., 2008; Sato & Shibasaki, 2013). Islands can also become incorporated in the adjacent floodplains as “docked islands” (Papacostas et al., 2008) while remaining hotspots of high arsenic groundwater. Figure 5 illustrates the occurrence of high As groundwater on Tay island on the Mekong floodplain (10°37′37″N, 105°22′27″E) with the majority of the wells containing more than 6.7 μmol/L As. The arsenic concentration is highest in the western part of the island, apparently reflecting the migration of the western channel in westerly direction. In addition to the high content of reactive organics and As‐containing Fe oxides of young sediments, hydraulic conditions support the high As concentration of groundwater on these islands. As pointed out by Sato and Shibasaki (2013), most groundwater discharge will occur at the sides of the river; there will be little lateral flow below the central part of the riverbed. The low vertical groundwater velocity allows ample time to build up a high As concentration from the kinetic release of As by reduction of Fe oxides.

Figure 5

Tay Island in the Mekong Delta, Vietnam, with high As groundwater. Based on Sato and Shibasaki (2013). Line indicates modeled cross section (Figure 6) with the position of the western channel at 2000 and 2500 years shown as crossing lines, the current channel position is present day at 3,000 years.

Figure 6

Westward channel migration of channel west of Tay Island. Panels (a)–(c): groundwater As after 2000, 2500, and 3000 years. Upper panel: groundwater flow pattern and RPVs flushed after 3,000 years.

A model was set up as a cross section through Tay Island, the Mekong, and the western channel. In the model, the western channel is moving westward in steps of 50 years, and for each step, recent sediment is added to its eastern margin. The resulting groundwater As concentration is shown after 2000, 2500, and 3000 years in Figure 6, and the corresponding channel position is indicated in Figure 5 with the present situation as the final stage. The recharge over the whole area is set to 0.15 m/year. The model aquifer thickness is set to 67 m with a downward flux of 0.015 m/year at the bottom. The model excludes the lower part of the 300 m deep aquifer described by Sato and Shibasaki (2013), to limit computation time. Initially, young sediment covers the area. After 2000 years of model time, elevated As concentrations have developed under the Mekong and the western channel. The velocity vectors (Figure 6) show the highest discharge to occur from the sides of the channel, and this water has relatively low As concentrations. There is upward flow around the river and channel that brings high As groundwater closer to the surface. The western channel then moves westward depositing recent sediment, and as shown for 2500 and 3000 years, the relatively near‐surface and high As groundwater also moves westward and persists over a considerable period of time (Figure 6), controlled by the adsorption capacity. This simulated high As groundwater corresponds to the observed high arsenic groundwater distributed parallel with the western channel (Figure 5). Overall, the association of islands and high arsenic groundwater is therefore the combined effect of the accumulation of recent reactive sediment and the discharge pattern bringing high As groundwater formed under the river channel upward.

Kandal Aquifer, Mekong Floodplain, Cambodia

The field area is situated between the Mekong and Bassac rivers about 20 km SE of Phnom Penh, Cambodia (11°31′28″N, 105°01′16″E), and has been studied extensively (Benner et al., 2008; Kocar et al., 2008, 2014; Polizzotto et al., 2008; Richards et al., 2017, 2018; Stuckey et al., 2015). The area is covered by wetlands, abandoned oxbows, and levees on top of a sedimentary sequence consisting of >50 m sandy aquifer overlain by 3–20 m red and gray clay deposits (Benner et al., 2008; Kocar et al., 2008). The aquifer sediments have a burial age of 7000–8400 years, but much younger (<2000 years) sediments have accumulated in oxbow lakes and ponds (Kocar et al., 2008).

Groundwater arsenic concentrations are in the range of 0.2–17 μmol/L (Figure 7) and are particularly high below the ponds that accumulate young sediment under anoxic conditions and in a downstream plume extending toward the Mekong river (Figure 7). Stuckey et al. (2015) showed in a series of sediment incubations that these permanent anoxic pond sediments readily released As, while the sandy aquifer sediment, or seasonal wetland sediment, released no arsenic.

Figure 7

Cross section through the Mekong floodplain, Kandal aquifer. Cambodia. Numbers in [ ] are μmol/L and in ( ) μg/L of As for the contoured intervals. Numbers in the small font indicate measured groundwater As concentration in μg/L. Based on Benner et al. (2008).

Using flow modeling, Benner et al. (2008) have estimated a vertical groundwater velocity of 0.04 to 0.4 m/year. The permeability distribution in the flow model was based on slug tests carried out in the clay layer. Subsequent excavations at the site have revealed fractures in the clay layer, which may cause a higher permeability than given by the slug tests (S. Benner, personal communication, October 24, 2017). This study did not include 3H/3He groundwater dating to directly measure flow velocities. Flow modeling suggests total groundwater travel times in the aquifer of 400–900 years (Benner et al., 2008) and that some 3–30 pore volumes have been flushed through the aquifer since its deposition (Polizzotto et al., 2008). However, Richards et al. (2017) recently measured a downward vertical groundwater velocity of 0.5 m/year, by using 3H/3He dating at Kien Svay, a few kilometers further SE. The groundwater dating by 3H/3He shows the groundwater to be much younger than predicted by the flow model of Benner et al. (2008). As a result, considerably more than 30 pore volumes must have recharged the aquifer since burial 7000–8400 years ago.

A 2‐D reactive transport model for this transect has previously been published by Kocar et al. (2014). In their model, organic matter degradation and Fe oxide reduction releasing As only occur in the young pond sediments, while in the sandy aquifer, sorption is the sole process affecting As. The model illustrates how high As groundwater emitted from the pond sediments is transported toward the river. The model of Kocar et al. (2014) uses the transport parameters given by Benner et al. (2008), implying that only 3–30 RPVs have gone through the aquifer. According to the model of Postma et al. (2016), this number of RPVs is not sufficient to flush arsenic from the aquifer. Sø, Postma, Vi, Pham, Kazmierczak, et al. (2018) pointed out that during flushing of the first 100 RPVs, high arsenic groundwater will persist in the aquifer. Radloff et al. (2017) applied their flushing model on this field site using the vertical groundwater velocity of 0.05 m/year of Benner et al. (2008). They found that it required 5000 years, close to the burial age of the aquifer, to leach As out of the aquifer to about 20 m of the total depth of 60 m. Both observations indicate an inconsistency between the low As concentrations observed in most of the aquifer (Figure 7) and the low estimated vertical groundwater velocity of 0.05 m/year.

We have constructed an alternative model by using the geochemical model of Postma et al. (2016) in a 2‐D model using the higher vertical groundwater velocity measured by using 3H/3He dating (Richards et al., 2017). The model is similar to that for the H‐transect in Vietnam, except that it extends to a depth of 60 m. The model was run for 7000 years, corresponding to the time from the initial deposition to the present burial age, and thereafter, a pond containing recent sediment was placed on the surface, and the model was run for an additional 700 years. The results (Figure 8) shows, in contrast to the model of Kocar et al. (2014), that even after 7000 years, high As concentrations are present in the deeper part of the aquifer, which corresponds to what was measured in the field (Figure 7). From the pond, a plume of As extends toward the river, similar to field observations (Figure 7) and the model of Kocar et al. (2014). Although the pattern of As distribution in the model (Figure 8) nicely fits the field observations (Figure 7), the concentration level in the model results is lower. This lower concentration is not unexpected since we have used the rate constants and geochemical sediment properties unaltered from the Red River site.

Figure 8

Conceptual model of the Kandal aquifer, Cambodia. Water has recharged for 7000 years at 0.175 m/year, followed by 700 years of recharge partially through recent pond sediments.

Discussion

Reactive Transport Modeling

We here present the first comprehensive 2‐D reactive transport model for evaluating groundwater arsenic concentrations in the floodplains of SE Asia. The same model was used to predict the behavior of arsenic at field sites thousands of kilometers apart, covering aquifer ages from 0 to 8000 years in different hydrogeological settings. It is based on the 1‐D geochemical model of Postma et al. (2016) developed for and calibrated on the data of the Red River floodplain presented by Postma et al. (2012). For transparency, we have refrained from adjusting geochemical rate parameters and sediment properties to the different field sites even though that would improve the model fits at each individual site. Our primary objective was to test the overall geochemical reaction network controlling the groundwater arsenic concentrations under a wide range of hydrogeological and sedimentary conditions. The satisfactory performance of the model under the variable conditions indicates that our model provides a good overall mechanistic description of the geochemical and hydrological processes controlling the groundwater arsenic concentration.

Combined Geochemical and Hydrogeological Controls

The main coupling between geochemical and hydrological processes is that for a given rate of arsenic release, the resulting groundwater As concentration depends on the velocity of the groundwater and the transport time, which are combined in the groundwater residence time. For a low groundwater velocity, the As release rate, in moles per unit time, is less diluted and the As concentration is higher, while for a fast groundwater velocity, the As concentration is lower. Due to the rate control on As release, the As concentration also becomes higher with increasing groundwater age (Richards et al., 2017; Stute et al., 2007) if this is the only process controlling the As concentration. Finally, as the most reactive organic carbon and Fe oxides are consumed over time, the rate of arsenic release decreases as time passes, while more RPVs pass through the sediment.

Sorption of As is an equilibrium reaction, and the chemical process itself is therefore not time dependent. Indirectly, the removal of sorbed arsenic is related to time via the number of RPVs that has passed, which is determined by the recharge rate. Furthermore, As sorption is related to time through the kinetic release of As by reductive dissolution of As‐containing FeOOH. The additional supply of released As delays the leaching of As from the aquifer, while the amount of FeOOH surface available for As sorption decreases. Thus, the two processes have the opposite effect on As leaching. The amount of FeOOH available for sorption and the site density of surface sites on the Fe oxide surface are accordingly important parameters for the rate of flushing arsenic out of aquifers.

To illustrate the combined effects of reductive dissolution of As bearing Fe oxides and sorption on the evolution of the groundwater As concentration, we have calculated their individual effects for the H‐transect example (Figure 9). Because the H‐transect model shows no variation in concentration along the x axis, data were extracted from the 2‐D model to show depth profiles for different model conditions. The result of reductive dissolution of As‐bearing Fe oxide without sorption is displayed in Figure 9a. Because organic carbon and Fe oxides are initially most reactive, the groundwater As concentration is highest after 100 years and shows the steepest concentration gradient over depth. The groundwater arsenic concentration and the gradient decrease over time with diminishing reactivity of the organic carbon and Fe oxide. The results of the full model, including both As release by reductive dissolution of FeOOH and As sorption to the diminishing FeOOH surface, are shown in Figure 9b and correspond to the 2‐D profile shown in Figure 2. The amounts of As(III) adsorbed or desorbed in the full model at different times are displayed in Figure 9c. At 100 years, most of the As(III) released by reductive dissolution becomes adsorbed and dissolved As(III) remains low. As the adsorbent fills with As(III), the aqueous As(III) concentration becomes higher. The As concentration increases over depth until it reaches the point where more As(III) can become adsorbed. Because the adsorbed As(III) concentration, expressed per liter of contacting groundwater, is much higher than dissolved As(III) (Postma et al., 2016), the aqueous As(III) concentration then becomes almost constant. In the upper part of the profiles, As(III) also starts to desorb (Figure 9c) as the kinetically controlled release rates become lower over time, and desorption of As(III) then contributes to the aqueous As(III) concentration. The interval with consistently elevated As(III) concentration (Figure 9b) passes through a maximum after 1200 years (Postma et al., 2016). Following this, the As(III) concentration decreases due to diminishing release rates but is still delayed by desorption of As(III). The amount of adsorbed As(III) is so high that it initially delays the increase in the groundwater As(III) concentration and later delays the decrease in As(III) concentration. The dynamics of sorption are further illustrated in Figure 9d, which considers an initial exchanger brought in equilibrium with a 10 μmol/L As(III) solution, and subsequently leached by As free water, corresponding to the model of Radloff et al. (2017). Traditional leaching fronts develop, but notably, their downward migration is slow, corresponding to a retardation of around 50. The model demonstrates that adsorption and desorption have a major influence on groundwater As(III) distribution over long timescales.

Figure 9

Effects of As release by reductive FeOOH dissolution and adsorption/desorption on groundwater As contents. Depth profiles extracted from the 2‐D model for the H‐transect, run at different times indicated by numbers indicating years. (a) Model with sorption removed. (b) Full model corresponding to Figure 2. (c) Desorption and adsorption from sediment expressed per liter of contacting groundwater. (d) Only sorption with initially 10 μmol/L As(III) in contact with sediment.

In summary, the release of As through reductive dissolution of FeOOH and subsequent sorption and desorption of As affect the groundwater As concentration in concert. The nature of the processes differs, giving a very different relationship to hydraulic conditions and evolution over time. Their combined effects have been quantified by using the model of Postma et al. (2016), which shows that the groundwater As distribution is a function of the number of RPVs flushed through the deposit since burial (Sø, Postma, Vi, Pham, Kazmierczak, et al., 2018). The model predictions agree with results from aquifers with different burial age and groundwater residence times and show that the highest As concentrations are found during the first 100 RPVs flushed through the aquifer, while about 200 RPVs are needed to bring groundwater arsenic down to the WHO level (Sø, Postma, Vi, Pham, Kazmierczak, et al., 2018).

A comparison of the H‐transect on the Red River floodplain, Vietnam, and the Araihazar aquifer, Ganges‐Brahmaputra floodplain, Bangladesh, shows two shallow sandy aquifers of almost similar burial age, respectively, 530 and 600 years. Both aquifers develop high groundwater As concentrations at shallow depth. The difference is that while the H‐transect is covered by a continuous fine‐grained top layer and displays a uniform As increase over depth (Figure 1), the Araihazar aquifer features sandy windows in the fine‐grained top layer with low As concentrations downstream from the windows (Figure 3). The difference in recharge rate between the two settings is about one order of magnitude (Figure 4), and the effect of the recharge rate on the As concentration illustrates the kinetic rate control on the release of As from the sediment. Accordingly, low groundwater arsenic can be present in young deposits as long as the recharge rates remain high. In terms of pore volumes, the sediment below the fine‐grained layers has typically been flushed by less than 100 RPVs, while the sediment below and just downstream the sand windows has been flushed by more than 300 RPVs (Figure 4).

Fine‐grained layers, which mostly reflect overbank flooding, on top of sandy aquifers, are a common feature of the floodplains of SE Asia. These layers are semiconfining and are often fractured due to desiccation in the dry season. 3H/3He dating of groundwater in sandy aquifers below such fine‐grained top layers shows slow downward vertical groundwater velocities. In the Araihazar aquifer, Bangladesh, Stute et al. (2007) reports vertical groundwater velocities down to 0.44 m/year, and Radloff et al. (2017) give velocities in the range 0.25–0.5 m/year, which with a porosity of 0.25 correspond to recharge rates of 0.06–0.125 m/year. Likewise, Postma et al. (2012) found vertical groundwater velocities of 0.25–0.5 m/year in sandy aquifers below fine‐grained top layers on the Red River floodplain. Below a fine‐grained cover (some 10 m of silt) in Holocene sands of W Bengal, McArthur et al. (2010) report vertical groundwater velocity of 0.5 m/year based on 3H/3He dating associated with As concentration up to 500 μg/L. Finally, Richards et al. (2017) measured a 3H/3He‐based vertical groundwater velocity of 0.5 m/year in an aquifer below a fine‐grained top layer on the Mekong floodplain, Cambodia. Nath et al. (2010) previously suggested that fine‐grained low permeability top layers on sandy aquifers may cause the development of high arsenic groundwater in young shallow aquifers, which is confirmed by the examples given here. A likely mechanism is that a low vertical groundwater flow velocity, in the range 0.25–0.5 m/year, allows time for the kinetics of organic matter degradation, to first develop anoxic conditions and then release arsenic through reductive dissolution of Fe oxides, building up a high arsenic concentration at shallow depth.

A massive clay layer, without fractures, on top of a sandy aquifer would not generate groundwater with As concentrations increasing over depth in the aquifer, because very little water is transported through such a low permeability clay layer. High dissolved arsenic concentrations could accumulate in the porewater of the clay layer, but diffusion of As out of the bottom of the clay layer into the aquifer would be the main transport mechanism. Diffusion would generate decreasing As concentrations with depth from the bottom of the clay layer, down through the sandy aquifer as the passing groundwater removes the arsenic. Generally, groundwater As concentrations increase with depth in sandy aquifers, suggesting that As diffusion from overlying fine‐grained layers is not an important process.

In shallow aquifers, fine‐grained confining layers may have an important role in limiting the number of RPVs that pass though the aquifer and, in combination with young burial age, promote the development of elevated As concentrations. In older Holocene aquifers, sediments with slower As release rates but with longer groundwater residence times may also generate high As groundwater. Thus, Sø, Postma, Vi, Pham, Kazmierczak, et al. (2018) reported an aquifer at 35 m depth with a burial age of 8000 years (Phuong Town) that developed high arsenic (2–4 μmol/L) groundwater due to a low vertical groundwater velocity resulting in less than 100 RPVs having flushed through the aquifer.

An additional example of hydrogeological control on the number of RPVs flushed and thereby the arsenic concentration is the flow situation below major river channels as elaborated upon in section 3.4. Here the rate of migration of the meandering river plays a role, and humans influence the process by building dykes along the rivers, which along the Red River has taken place for at least 1000 years. Dyke building confines the river meandering and may thereby result in building up high groundwater arsenic because of accumulating young reactive sediment between the dykes and increased groundwater residence times under the river. The arsenic concentrations will depend on the vertical flow rates. Thus, Stahl et al. (2016) found a low groundwater recharge rate in the river bed of the Red River at Van Phuc that resulted in a high groundwater As concentration. In contrast, 6 km upstream the Red River at Nam Du, where pumping for Hanoi's water supply causes downward groundwater velocities of up to 19 m/year, the same sediment with the same high As release rate yields low As groundwater because of the large dilution of the released As (Postma et al., 2017).

Finally, the Cambodia example shows how sediments with a different burial age located adjacently may generate very different groundwater arsenic concentrations. It indicates that the geological structure of the floodplains, with emphasis on intersecting channels that place sediments of various ages side by side, is a likely cause of the variability of the groundwater As content (McArthur et al., 2010; Weinman et al., 2008). The use of remote sensing might be a helpful tool to understand the floodplain structure and identify the youngest sediments where the highest As concentrations can be expected (Papacostas et al., 2008; Quicksall et al., 2008). The Cambodia example also shows how groundwater from young sediment generating high As can be transported into older sediments; this complexity, generated when layers of different age are stacked into a sequence, will be elaborated in the next section.

Geological Evolution and Groundwater Arsenic

So far we have demonstrated how present‐day variability in geochemical and hydrogeological parameters influences the As concentration of groundwater. This discussion, however, neglects the changes that occur over geological time. Sedimentary flood plain basins subside and sediments accumulate or erode, river channels shift position, and the hydrogeological conditions change as well. In a first attempt to analyze how the geological evolution in a flood plain basin may affect the groundwater arsenic concentration, we have analyzed a theoretical scenario of a floodplain with shifting channels in an accreting sedimentary sequence. The scenario schematically follows the hydrogeological evolution on the Red River floodplain upstream Hanoi over the last 6000 years (Kazmierczak et al., 2016).

The scenario modeled comprises three time steps: The first step (0–2400 years) considers a Holocene floodplain with a centrally located main channel placed in sediment with a relatively high hydraulic conductivity (Figure 10). Initially the entire model is 16.5 m thick, the lower 4 m are Pleistocene sand and clay, making the Holocene deposits 12.5 m thick, and the model transect is 17 km across. The Holocene deposits are hydraulically separated from the underlying Pleistocene sand by the clay layer, but there are some connections between the two aquifers. The recharge is constant, and the flow is controlled by the head defined for the rivers and channels defined as drains, as well as the constant head boundary at the bottom of the model. Using this initial hydrogeological configuration, the reactive transport model is run for 2400 years with the steady‐state flow pattern shown in Figure 10a and the resulting As distribution in the groundwater displayed in Figure 11a. The model result shows groundwater flowing to the main channel as well as to some minor streams (Figure 10a). Generally, the As concentration increases over depth (Figure 11a) but is particularly high beneath the main channels, due to low groundwater flow, as discussed in the previous Tay Island section. In the second step (5400–2400 years), the upper 2.75 m of the sequence has become eroded while 4.75 m of floodplain sediment are added on top increasing the total thickness to 18 m. Meanwhile, the channel shifts location to the left and erodes into the underlying deposits (Figure 10b). For this new hydrogeological situation, with new defined drains and associated heads, the reactive transport model is now run for 5400–2400 = 3000 years, starting with the geochemical properties in the sediments surviving erosion while recent sediment is added on top. During these 3000 years, a new flow pattern exists (Figure 10b) and the groundwater arsenic distribution is modified as well (Figure 11b), with high As below the new position of the channel and remnants of high As groundwater below the previous channel position. The third step (5400–6000 years; Figure 10c & Figure 11c) brings us up to the present day situation. An additional 4 m of sediment is added on top of the floodplain, to a total sequence of 22 m, while the main channel now shifts to the right‐hand side of the profile and cuts into the underlying deposits. Again, the heads of the drains are changed leading to a new flowfield. And again, the geochemical properties of model cells that are not affected by erosion are retained and defined the initial state of the subsequent reactive transport model, with recent sediment added on top before the start of the last 600 years of reactive transport modeling. The resulting groundwater arsenic distribution shows particularly high concentrations in the youngest channel sediments. Comparison of the arsenic distribution at 6000 years with the number of RPVs flushed since burial (Figure 11d) shows a good overall correspondence in the upper layers while at medium depth the arsenic concentrations are higher than expected from the RPV distribution probably due to downward transport of high arsenic water into older sediments. Overall, this modeling exercise, though simplified compared to the actual system, illustrates the complexity in As distribution generated when a sequence of river sediments accumulates.

Figure 10

The geological development in the three‐step model. Lower three panels: Sedimentation and channel migration as reflected by hydraulic conductivity, and groundwater flow induced by the recharge and the drains marked by green lines on the curves showing the position of the water table. Upper panel: Sediment age distribution at 6000 years.

Figure 11

Modeled As distribution after 2400, 5400, and 6000 years as indicated in Figure 10. Upper panel shows the RPVs since burial.

Conclusions

Neither hydrogeological nor geochemical processes by themselves can explain the arsenic distributions commonly observed in aquifers on the SE Asian floodplains. It is the interplay between the burial age of the deposits, the kinetic control of arsenic release from sediment, the sorption/desorption, the groundwater residence times, and the recharge rate that determine both the release of arsenic from the sediment and the flushing of arsenic from the aquifer. These interactions can only be quantified by using reactive transport modeling.

We have shown that our 2‐D reactive transport models using the geochemical model of Postma et al. (2016) can reproduce the groundwater arsenic distribution in four well‐studied recent settings throughout SE Asia in aquifers with highly diverse hydrogeological conditions and a wide range in aquifer sediment burial age.

Using our reactive transport model for modeling the effect of floodplain evolution over 6000 years on the groundwater arsenic content indicates overlapping timescales between the rates of channel migration and the dynamics of arsenic mobilization. The result is evolving complex patterns of groundwater arsenic distribution, with patches of high arsenic groundwater as chemical afterimages of vanished rivers. Features that cannot be explained by the observable present state of the system.

A quantitative understanding of the spatial distribution of arsenic in deeper groundwater requires the understanding of the hydrogeological evolution of the stacked geological sequences, which imposes a major challenge both in terms of the collection of the needed field data and reactive transport modeling capacity.