Journal of Hazardous Materials 280 (2014) 685–695

Contents lists available at ScienceDirect

Journal of Hazardous Materials journal homepage: www.elsevier.com/locate/jhazmat

Reactive transport modeling of 90 Sr sorption in reactive sandpacks Jun Yin a,b,∗ , Sung-Wook Jeen c,d , David R. Lee e , K. Ulrich Mayer a a

Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, BC, Canada Ministry of Forests, Lands and Natural Resource Operations, Prince George, BC, Canada c Department of Earth and Environmental Sciences, Chonbuk National University, Jeonju, Jeollabuk-do 561-756, Republic of Korea d The Earth and Environmental Science System Research Center, Chonbuk National University, Jeonju, Jeollabuk-do 561-756, Republic of Korea e Atomic Energy of Canada Limited, Chalk River Laboratories, Chalk River, ON, Canada b

h i g h l i g h t s • • • •

Evaluation of reactive sandpacks for in situ remediation of radioactive 90 Sr. Quantitative analysis of in situ column experiments with fast flow rates. Results show that velocity-dependent kinetic sorption is important. Sensitivity analysis for the design of reactive sandpacks in dewatering projects.

a r t i c l e

i n f o

Article history: Received 22 February 2014 Received in revised form 9 July 2014 Accepted 28 July 2014 Available online 26 August 2014 Keywords: Strontium-90 (90 Sr) Reactive transport modeling Kinetic sorption Dewatering well

a b s t r a c t Strontium-90 (90 Sr) is one of the most problematic radioactive contaminants in groundwater at nuclear sites. Although 90 Sr is retarded relative to groundwater flow, it is sufficiently mobile and long-lived to require treatment in many hydrogeological settings. A detailed study was performed on the practicality of using granular clinoptilolite as a sandpack around groundwater wells where groundwater is contaminated with 90 Sr and the water table must be lowered. The effectiveness of the reactive sandpack concept and the mechanisms controlling 90 Sr attenuation was investigated by numerical analysis of data obtained from four in situ column experiments. The experiments spanned the range of pore-water velocities that would occur during radial flow through granular clinoptilolite sandpacks. A kinetic sorption model was required to adequately reproduce the experimentally observed 90 Sr behavior. Calibrated first-order kinetic rates were correlated with pore-water velocities. After calibration, three sorption models were used to simulate 90 Sr attenuation for four hypothetical pumping scenarios. Results show that a velocity-dependent kinetic model accurately simulates the observed early breakthrough for high pore-water velocities. The results indicate (1) that reactive sandpacks have good potential for in situ remediation and construction dewatering and (2) that quantitative modeling can aid in the design and application of this novel technique. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Strontium-90 (90 Sr) is one of the problematic radioactive contaminants in groundwater at many nuclear sites. The half-life of 90 Sr is 29.1 years and the daughter isotope 90 Y will decay to 90 Zr with a half-life of 64.1 h. For example, at Chalk River Laboratories, Chalk River, Ontario, Canada, 90 Sr was released into the

∗ Corresponding author at: Ministry of Forests, Lands and Natural Resource Operations Water Stewardship Division 1044-5th Ave Prince George, BC Canada V2L 5G4. Tel.: +1 250 565 6440; fax: +1 250 565 6629. E-mail addresses: [email protected], [email protected] (J. Yin). http://dx.doi.org/10.1016/j.jhazmat.2014.07.073 0304-3894/© 2014 Elsevier B.V. All rights reserved.

local aquifer during the 1950s [1]. Although 90 Sr migration was naturally retarded by sorption onto mineral surfaces, the groundwater plume eventually reached a location in the groundwater flow system where it was about to impact the biosphere [2]. The geochemical mechanisms of 90 Sr attenuation in natural aquifers are reasonably well understood [1,3,4]. Subsurface retention of 90 Sr can be vastly improved with the use of specific sorptive materials [5,6]. The potential of the zeolite clinoptilolite to retard the migration of 90 Sr in groundwater was proposed by Lee and Hartwig [7]. In 1998 a permeable reactive barrier (PRB), using granular clinoptilolite as the reactive medium, was installed at Chalk River Laboratories (CRL) to remediate 90 Sr contaminated groundwater [7]. As the reactive constituent of the PRB, clinoptilolite has proven

686

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

to be an effective ion-exchange reactant and has demonstrated that the 90 Sr contaminated groundwater can be treated in situ under natural groundwater flow conditions [8]. This paper focuses on another application of clinoptilolite, using it as a reactive sandpack for dewatering wells in contaminated ground. A reactive sand pack would be an “active” remediation method, as opposed to the “passive” PRB systems in which groundwater migrates through the treatment material under natural hydraulic gradients. Sandpacks are routinely used in construction dewatering [9] and the concept of reactive sandpacks may be suited for dewatering projects where the groundwater is contaminated. If clinoptilolite can be used to retain subsurface contamination in situ around the screens of dewatering wells, the cost for the storage and treatment of the pumped groundwater can be substantially reduced. In situ reactive sandpacks might also be useful as preliminary treatment for above-ground treatment systems. The use of granular clinoptilolite as a reactive sandpack material for dewatering well screens was first applied during the construction dewatering around the sheet pilings used to build the PRB at CRL. In this case, the predominantly native sandpack was replaced by an appropriately sized granular clinoptilolite sandpack. The physical efficiency of the dewatering wells was unaffected, but the actual effectiveness of the clinoptilolite to decontaminate the pumped water was not quantified [11]. More recently, a laboratory study of reactive sandpacks was conducted [10,12] and the results were promising. One of the major challenges of using the reactive sandpack approach is that reaction kinetics may reduce remediation efficiency as the pore-water velocity increases near the well screen. Previous investigations showed that the average linear groundwater velocity at the Chalk River PRB site is 0.2 m d−1 [8]. This value is several orders of magnitude lower than pore water velocities near pumping wells. While the concept of reactive sandpacks appears useful and the sorption capacity of clinoptilolite is considerable, it is unclear if reaction kinetics will be rapid enough for satisfactory removal of contaminants. The purposes of this study are to evaluate kinetic limitations on reactive sandpacks as a means of remediation for 90 Sr contaminated groundwater and to develop a reactive transport modeling tool to interpret experimental data and to aid in the design of field applications. In this paper, the results of four in situ column experiments on 90 Sr sorption for different flow velocities are presented, two sorption models are calibrated against the observed 90 Sr breakthrough curves, the calibrated models are used to simulate 90 Sr attenuation in reactive sandpacks of hypothetical dewatering wells and, finally, simulation cases are tested to analyze the sensitivity of 90 Sr breakthrough as a function of pumping rate and sandpack thickness. 2. In situ column experiments In-well, field-column experiments were conducted to test the attenuation capacity of reactive sandpacks for a range of flow

90 Sr

velocities. The experiments were conducted in the shallow portion of the 90 Sr plume at WMA “A” in the Chalk River Laboratories to maintain natural geochemical conditions and to avoid handling contaminated water in the laboratory. Four 10-cm long columns were set up and tested to evaluate the effectiveness of reactive sandpacks. More details of the column setting and sample analysis can be found in the supplementary data. Before deploying the columns in the field, tracer (600 mg L−1 NaCl) tests were conducted in the laboratory using similar flow rates as used for the in situ field experiments. Electrical conductivity (EC) of the effluent was measured and recorded. After completion of the tracer tests, deionized water was used to flush the columns until the EC returned to the pre-tracer test levels. The physical transport parameters for the material, specifically porosity and longitudinal dispersivity, were calibrated using the modeling package CXTFIT [13]. 3. Experimental results and model calibrations 3.1. Column experiments 3.1.1. Laboratory Tracer tests using NaCl were conducted in the laboratory before exposing the columns in the field to groundwater containing 90 Sr. Table 1 shows the physical properties of the four columns and the calibrated transport parameters. The four observed and calibrated EC curves are plotted in Fig. 1(a)–(d). Overall, no significant ratelimited mass transfer was observed. Transport could be simulated using the classic advection-dispersion equation without dual or multi-domain effects. Calibrated porosities (0.40–0.46) were similar to the values calculated by measuring column mass before and after saturation with water (0.43–0.44). Longitudinal dispersivities varied slightly within a range of 0.95–1.45 cm. The data proved that the columns were performing as intended (i.e., no evidence of turbulent flow or channeling). 3.1.2. Field Four field columns were installed below ground in a well that was screened in groundwater containing 90 Sr. This was done to expose the reactive material to the same chemical conditions as would occur if a reactive sandpack were used at this location. During field operation, column flow rates and inflow 90 Sr concentrations fluctuated as shown in Fig. 2(a)–(d). Major ion chemistry also varied slightly between the columns and over time (Table S1, Supplementary Data). Flow rates were likely affected by the precipitation of iron oxides and possibly gas exsolution in the column and sampling lines. For example, the sharp decrease in flow rate after day 31 in column C4 (Fig. 2d) was likely due to clogging of the column by precipitation of iron oxides. The formation of iron oxides was supported by the observed reddish color in the treatment material after sectioning the column. These results indicate that oxygen was entering the columns during the experiment or that colloidal iron was already present in the groundwater. The

Table 1 Physical properties of the four in situ columns when tested in the lab using NaCl as a conservative tracer. Column

C1

C2

C3

C4

Pumping rate (pore water velocity)a Operation period (days) Dry mass (g) Bulk density (g cm−3 ) Pore volume (cm3 ) Porosity (gravimetric) Porosity (tracer test) Dispersivity (cm) (tracer test)

17 mL min−1 (28 m d−1 ) 49 175.65 0.89 84.25 0.43 0.40 0.95

28 mL min−1 (48 m d−1 ) 54 175.76 0.89 84.45 0.43 0.40 1.22

48 mL min−1 (81 m d−1 ) 55 175.65 0.89 85.83 0.44 0.46 0.98

94 mL min−1 (155 m d−1 ) 50 (31)b 173.67 0.88 87.17 0.44 0.45 1.45

a b

In the field, the actual velocities in the experiments are different than the designed flow velocities. When this column was used in the field the velocity dropped to nearly zero in Column C4 after this time period.

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

b) Column C2 (at 20 mL min-1)

1.0

1.0

0.8

0.8

0.6

0.6

EC (C/C 0)

EC (C/C 0)

a) Column C1 (at 21 mL min-1)

0.4

0.2

0.0

200

300

400

500

600

0.4

0.2

Measured Fitted 100

0.0

700

100

200

Time (seconds)

400

500

600

700

d) Column C4 (at 114 mL min-1)

1.0

1.0

0.8

0.8

0.6

0.6

EC (C/C 0 )

EC (C/C 0 )

300

Time (seconds)

c) Column C3 (at 132 mL min-1)

0.4

0.2

0.0

687

0.4

0.2

50

100

150

200

0.0

250

50

TIme (seconds)

100

150

200

250

Time (seconds)

Fig. 1. Measured and simulated electrical conductivities (EC) for NaCl tracer tests in four columns.

average column flow rates, based upon the accumulated effluent volume, were 17, 28, 48, and 94 mL min−1 , respectively. This corresponds to average linear groundwater velocities of 28 (Fig. 2a), 48 (Fig. 2b), 81 (Fig. 2c), and 155 m d−1 (Fig. 2d). The reason for the fluctuating influent 90 Sr concentration (Fig. 2(e) and (f)) appeared to be due to natural variations in groundwater as described above, possibly aggravated by transient changes in groundwater flow directions. Despite the strong sorption capacity of the clinoptilolite, breakthrough was observed in all four columns within the duration of the experiments. Fig. 3 displays field observations and the results of modeling with a kinetic and an equilibrium model (described below). Clearly in column C4 the field observations were unexplainable with effluent values rising to nearly 2 × 10−12 mol L−1 and then falling erratically to much lower values. Although scattered, breakthrough of 90 Sr was correlated to the flow rates of the columns, indicating that the sorption was rate-limited. In columns C3 and C4, the fluctuations of the effluent 90 Sr concentrations were correlated to the 90 Sr concentration changes in the influent. This provided further evidence for kinetic limitations at the faster flow rates. Fig. 4a shows the inventory of 90 Sr along the lengths of the four field columns with the inventory correlated to the mass loading. Fig. 4(b)–(e) show the match between the measured and simulated concentrations in the clinoptilolite. The relatively even distribution

of 90 Sr along the length of columns C3 and C4 suggests that uptake took place throughout the column (Fig. 4a). The lower adsorbed 90 Sr mass near the inlet of column C3 may be attributed to the sampling variability. In columns C1 and C2, the adsorbed 90 Sr concentrations were concentrated near the inflow end and decreased along the column (Fig. 4a). Both the lower adsorbed mass in Fig. 4a and the much lower breakthrough concentrations (Fig. 3(a) and (b) at 28 and 48 m d−1 , respectively) show that kinetic effects played less important roles at the two lower flow rates. In general, the maximum sorption capacity had not been reached at the end of the experiment in all four columns as shown by lower inventories at the outflow ends.

Table 2 Mass balance calculations of 90 Sr (mol) in the four experimental columns used in the field to mimic various positions in a hypothetical reactive sandpack of granular clinoptilolite.

C1 C2 C3 C4 a

Influent

Effluent

Sorbed

% Errora

2.11E−09 3.85E−09 8.55E−09 9.50E−09

3.06E−11 3.53E−10 2.52E−09 3.53E−09

1.32E−09 3.20E−09 4.68E−09 6.15E−09

36.3 7.7 15.8 1.9

% Error = abs [(cout v + csorbed vb ) − cin v]/cin v × 100.

688

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

140

140

a) Column C1

120

Pumping rate (ml min-1)

Pumping rate (ml min-1)

120 100 80 60 40 20 0

0

10

20

30

60 40

0

0

10

20

Time (days)

140

c) Column C3

30

40

d) Column C4

120

Pumping rate (ml min-1 )

Pumping rate (ml min-1 )

80

40

Time (days)

120 100 80 60 40 20

100 80 60 40 20

0

4E-12

10

20

30

Time (days)

40

0

50

0

10

20

30

40

Time (days)

4E-12

50

f) Columns C3 & C4

Sr-90 concentrations (mol L-1)

Sr-90 concentrations (mol L-1)

e) Columns C1 & C2

3E-12

3E-12

2E-12

2E-12

1E-12

0

100

20

140

0

b) Column C2

1E-12

0

10

20

30

40

0

50

Time (days)

0

10

20

30

40

50

Time (days)

Fig. 2. Measured pumping rates and inflow 90 Sr gross beta in four columns placed in a well that had been screened in the groundwater plume.

The total 90 Sr mass recovery is shown in Table 2 and the error is calculated as:

 % errror =

  c V− cin V  out

csorb Vb +

cin V

(1)

where csorb is the sorbed 90 Sr concentration (mol g−1 ); V is the accumulated influent/effluent volume (mL) between two samplings; b is the bulk density (g cm−3 ); and cin and cout are the inflow and outflow 90 Sr concentrations (mol L−1 ), respectively. Calculated mass balance was generally good with errors ranging from 2 to 16% for

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695 -1

689

-1

a) 28 m d (Column C1)

b) 48 m d (Column C2)

1E-13

8E-14

Sr-90 concentrations (mol L-1 )

Sr-90 concentrations (mol L-1 )

Measured Equilibrium Kinetic

4E-13

3E-13

6E-14

2E-13

4E-14

1E-13

2E-14

0

0

10

20

30

0

40

0

10

20

Time (days)

-1

40

50

-1

c) 81 m d (Column C3)

d) 15 5 m d (Column C4)

Sr-90 concentrations (mol L-1 )

2E-12

Sr-90 concentrations (mol L-1 )

2E-12

1E-12

0

30

Time (days)

1E-12

0

10

20

30

40

50

0

0

10

Time (days)

20

30

40

50

Time (days)

Fig. 3. Measured 90 Sr concentrations and simulated 90 Sr breakthrough curves in four field columns.

columns C2, C3, and C4, while column C1 had a larger error of 36%. Given that column C1 has the slowest pumping rate and therefore the least mass input, the error is likely due to the measurement uncertainties associated with lower concentrations. 3.2. Model calibration To interpret processes controlling 90 Sr sorption and provide design parameters for reactive sandpacks, sorption models were calibrated against data from the column experiments. In this study, the multi-component reactive transport model MIN3P [14] was used to conduct the calibration and hypothetical field scale simulations. The model solves the reaction advection-dispersion equations using the global implicit solution approach [15] and has been used in a number of previous reactive transport studies [e.g., 14,16,17]. Since kinetic behavior was observed in all four columns, rate-limited sorption was implemented for the present study. A first-order kinetic sorption equation was used in the model to express the source and sink term in the advection-dispersion equation:

∂Csorb = ˛(S − Csorb ) ∂t

(2)

where Csorb is the sorbed concentration (mol g−1 ); ˛ is the firstorder kinetic rate coefficient (d−1 ); and S is the sorption capacity of the chemical component (mol g−1 ), which is the maximum sorbed concentration when equilibrium is reached. The sorption capacity S can be expressed using a linear sorption model or non-linear Freundlich or Langmuir models. Model calibrations using all three sorption models were conducted, and the differences were minor (data not shown), suggesting that a more complex nonlinear modeling approach is not warranted. As a result, only the linear sorption model was used for the analysis of the experimental data. The sorbed 90 Sr was normalized to the pore water volume in the simulation using  Cˆ sorb = (1 − )s Csorb

(3)

where  Cˆ sorb is the sorbed concentration normalized to the pore water (mol L−1 ); s is the solid density of the material (g cm−3 ); and  is the porosity (−). Substituting Eq. (3) into (2) and using the linear sorption model (Kd = S/Caq , where Kd represents the distribution coefficient (mL g−1 ) and Caq represents the aqueous concentration (mol L−1 )) yields:

∂Cˆ sorb =˛ ∂t



Kd (1 − ) s Caq − Cˆ sorb 



(4)

690

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

a) Columns C1, C2, C3 & C4 C1 C2 C3 C4

Sr-90 concentrations (mol L-1 )

8E-08

6E-08

4E-08

2E-08

0

2

4

6

8

Distance from the inlet (cm)

1.5E-07

1.5E-07

1E-07

1E-07

Measured Equilibrium Kinetic

5E-08

0

c) 48 m d-1 (Column C2)

Sr-90 concentrations (mol L-1 )

Sr-90 concentrations (mol L-1 )

b) 28 m d-1 (Column C1)

5E-08

0

2

4

6

8

0

10

0

Distance from the inlet (cm)

1.5E-07

4

6

8

10

Distance from the inlet (cm)

1.5E-07

e) 155 m d-1 (Column C4)

Sr-90 concentrations (mol L-1 )

Sr-90 concentrations (mol L-1)

d) 81 m d-1 (Column C3)

1E-07

1E-07

5E-08

0

2

5E-08

0

2

4

6

8

10

0

0

Distance from the inlet (cm) Fig. 4. Measured and simulated 90 Sr concentrations for solid samples in four columns. The sorbed inventory determined in all 4 columns at the end of the in situ experiment.

2

4

6

8

10

Distance from the inlet (cm) 90

Sr was normalized to the pore water volume. Panel (a) shows the

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

3.3. Modeling results The calibrated breakthrough curves for the equilibrium model are shown in Fig. 3. As mentioned in Section 3.1, in column C4, the sharp decrease in flow velocity was likely caused by pore space clogging due to iron oxide precipitation. Although the decreased flow velocity was incorporated, clogging would also affect the characteristics of the sorption sites (surface area, site density, etc.) but this process was not investigated in the present study. Accordingly, the last eight measurement points from day 31 onwards (Fig. 3d) were given less weight in the parameter estimation. Using the equilibrium approach, the variations of 90 Sr influx were significantly buffered by sorption and desorption in all four columns. The contaminant was fully retarded at the beginning, and  breakthrough occurred at a later time. The four optimized K d values are 31,400 in column C1, 39,000 in column C2, 47,600 in column C3, and 45,600 in column C4. The corresponding Kd values are 15,200 mL g−1 , 18,900 mL g−1 , 23,500 mL g−1 , and 22,800 mL g−1 , respectively. These results can be compared to previous studies on Srsorption. Using similar sorption materials, Rabideau et al. [20] obtained 2965 mL g−1 in batch sorption isotherm experiments and 2045 mL g−1 in column experiments with slower pore water velocities. Fuhrmann et al. [21] obtained 630 mL g−1 and 685 mL g−1 in two column experiments also with slower pore water velocities. Samper et al. [22] calibrated 1114–1823 mL g−1 for Sr attenuation in compacted Ca-bentonite using inverse modeling. Comparing the distribution coefficients from these studies to the present results reveals that previously obtained Kd values were significantly lower. However, calcium levels in the studies by Rabideau et al. [20] and Fuhrmann et al. [21] were around 100 mg L−1 , substantially higher than the ambient Ca2+ concentration in the Chalk River site groundwater, which is about 10 mg L−1 (Table S1, Supporting Data). Priebe [10] conducted column experiments also using low calcium concentrations of 20 mg L−1 , and obtained Kd values comparable with our results (5000–32,000 mL g−1 ). The effect of competitive sorption involving Ca2+ in radionuclide transport has been well studied [e.g., 1,23,24]. With increasing Ca2+ concentration, sorption capacity for Sr2+ tends to decrease due to competitive sorption. The Kd values reported here must therefore be considered apparent Kd values valid only for the groundwater composition where the experiments were conducted. Ambient Ca2+ during the first set of experiments (C3 & C4, 7.5–8.2 mg L−1 ) was slightly lower than the

3.0 2.5 -1

Kineticrate(day )

Therefore, in the model calibration, the calibrated distribution   coefficient (K d ) becomes ((Kd (1 − ))/)s . Note that K d is dimensionless. For the simulations, the one-dimensional column was discretized with a spatial resolution of 1 cm. In order to account for transient boundary conditions (Fig. 2e and f), measured pumping rates and 90 Sr concentrations were applied at the column inlet using a specified flux (third-type) boundary condition. Outflow boundary conditions were specified hydraulic head for water flow and a free exit for solute transport. The model-independent parameter estimation program PEST [19] was used to optimize the sorption parameters for equilibrium and kinetic model implemen tations. In the calibration, the K d values were first estimated for the equilibrium model. When calibrating the kinetic models, the  K d values (estimated from the equilibrium model) were fixed and only the kinetic parameter ˛ was optimized. Considering the uncertainties associated with the column sectioning and solid concentration measurements, as well as the resulting mass balance errors (Table 2), the model calibrations were performed against the aqueous 90 Sr measurements. However, adsorbed 90 Sr concentration distributions were used to verify the calibration, keeping in mind the mass balance errors.

691

α =9.566v 2 R =0.793

2.0

-1.748

1.5 1.0 0.5

0

5

10

15

-1

Porewatervelocity(mday )

20

Fig. 5. Regression curve between the pore water velocity and kinetic rate.

Ca2+ in the second set of experiments (C1 and C2, 9.4–10.9 mg L−1 ) (Table S1, Supplementary Data). Therefore, the apparent Kd values in columns C3 and C4 were slightly higher than the values in columns C1 and C2, further corroborating the dependence of 90 Sr sorption on Ca2+ concentrations. Due to limited availability of major ion data; however, it is difficult to quantify the dependence of Kd values on Ca2+ concentrations. In columns C1 and C2, both simulations showed similar results between the kinetic approach and the equilibrium Kd approach, indicating near equilibrium behavior for this range of flow velocities. However, the observed 90 Sr concentration oscillations in column C2 were slightly better reproduced using the kinetic approach. Incorporating rate-limited sorption significantly improved the ability of the model to match the observed breakthrough curves for the faster flow velocities (columns C3 and C4). Using the kinetic model, early breakthrough could be reproduced for both columns. In addition, aqueous 90 Sr fluctuations at early stages were simulated, although with dampened amplitude (Fig. 3c and d). Simulated fluctuations diminished toward the end of the simulation, attributed to desorption of adsorbed 90 Sr during periods of low 90 Sr inflow concentrations. Optimized first-order kinetic rate constants in the four columns were 2.133 d−1 in column C1, 0.614 d−1 in column C2, 0.100 d−1 in column C3, and 0.135 d−1 in column C4. Dai et al. [25] used half sorption reaction times and Damköhler numbers to determine when a kinetic sorption model should be used to simulate solute transport. Using the same approach, the calculated half reaction times in the four columns range from 7.8 h to 123.2 h: significantly longer than the mean residence time of 90 Sr in the column. The calculated Damköhler numbers range from 7.6 × 10−3 to 8.7 × 10−5 , indicating that all the reaction rates are much slower than the transport rates (Damköhler number  1). Both criteria confirm that the use of a kinetic sorption model is required. Further examination reveals that the first-order rate constants are negatively correlated to pore water velocities and follow a power–law (Eq. (5) and Fig. 5): ˛ = 9.566v−1.748

(5)

where v is the pore water velocity (m d−1 ). On the microscale, sorption kinetics are a function of attachment and detachment processes at the mineral surface and are dependent on material characteristics including mineralogy and intra-grain pore distributions [26]. However, on the macro-scale, the magnitude of the sorption kinetics is dependent on pore water flow velocities, which has been observed in numerous previous

692

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

Table 3 Parameters of the four field scale simulation cases. Parameters

Case 1

Well radius (m) Aquifer thickness from the well screen (m) Aquifer depth (m) Initial 90 Sr concentration (mol L−1 ) 90 Sr concentrations in influent (mg L−1 ) Sand pack porosity Aquifer porosity Longitudinal dispersivity in both aquifer and sandpacks (cm) Reactive sand pack thickness (m) Pumping rate (m3 s−1 )

0.15 3

a b

Case 2

Case 3

Case 4

10 1.0 × 10−20 3.43 × 10−7 a 0.44 0.35 10

1

1 −2 b

1.5 × 10

2 −2

3.0 × 10

2 −2

1.5 × 10

3.0 × 10−2

This is equivalent to 875 Bq L−1 , which is the averaged 90 Sr concentration in the plume where the in situ column experiments were conducted at the Chalk River site. Typical pumping rate at the Chalk River site [10].

studies. For example, El-Kamash [18] studied the Sr2+ sorption on zeolite and observed decreased breakthrough times with increased flow rates. The uneven distribution of the influent inside the column due to the reduced contact time was used to explain the observations. An inverse correlation between the kinetic rate and the water velocity was also observed by Kim et al. [27]. Maraqa [28] reviewed several previous studies and concluded that the relationship between the mass transfer coefficient and pore water velocity follows a power–law relationship. The actual relationship, however, depends on the investigated media. In general for higher flow velocities, the residence time of pore water in the column becomes shorter, as does the time for diffusion of 90 Sr into intra-granular pores, and the 90 Sr is therefore more mobile. Observed and simulated adsorbed 90 Sr concentrations are directly compared in Fig. 4b to e. Similar to the aqueous 90 Sr concentrations, incorporating the kinetic approach improved the simulation results. The most significant differences between two models were in the columns with the higher flow velocities (C3 and C4). In column C3, total sorption mass was reduced by 15.5% using the kinetic sorption model. Lower 90 Sr concentrations near the column inlet were not reproduced in both models, although transient boundary conditions were considered. The reason for this behavior remains unknown, but is probably due to sampling variability or the heterogeneity created during the column packing. In column C4, the reduced sorption mass due to the kinetic effects was 19%. In addition to the mass balance, the 90 Sr distribution in the column was also improved significantly. The optimization of kinetic rates for the columns with low flow velocities (C1 and C2) showed very limited differences between the equilibrium model and the kinetic model, suggesting that kinetic effects are less important. 4. Evaluation of reactive sandpacks for 90 Sr treatment in the field In a dewatering well, pore-water velocities increase with decreasing distance from the well screen, resulting from the reduced cross sectional area. Based on the calibration results, sorption kinetics in the sandpack will be affected by velocity variations and the efficacy of treatment will be substantially reduced near the well screen. Therefore, when approaching the well screen, the slower sorption kinetics associated with the higher velocity will reduce the reaction rate and rapid breakthrough may occur. To evaluate the effect of kinetic sorption on the performance of reactive sandpacks under radial flow conditions, a sensitivity analysis was conducted for a hypothetical reactive sandpack field installation. Key parameters including thickness of the sandpack and pumping

rates were varied to formulate four simulation cases designed to analyze the sensitivity of 90 Sr treatment to these parameters. The simulations were conducted using a radial coordinate system [29,30] for a homogeneous model domain. It was assumed that pH and major cation concentrations in the groundwater are at steady state, and a constant Kd was used for all four cases. For simplicity, the average Kd from the four column experiments (Kd = 20,092 mL g−1 ) was used. The velocity-dependent kinetic sorption model, as calibrated from the in situ column experiments, was employed and the rate in each computational volume was calculated using the power–law relationship defined in Eq. (5). In addition to the velocity-dependent kinetic model, the constant rate kinetic and the equilibrium models were also used to evaluate the sensitivity of the results to the model formulation. For the constant rate model, the first order rate coefficient was set to 0.135 d−1 . This corresponds to the calibrated kinetic rate in column C4, where the pore water velocity was similar to the velocity used in the hypothetical test cases. All cases were run for a time period of 800 days. In all simulations, the dewatering well had a well radius of 0.15 m and the well was assumed to be screened over the entire depth of the aquifer (10 m). In Case 1, a pumping rate of 0.015 m3 s−1 was used [10] and the thickness (radius) of the sandpack installed around the screen was 1 m. By varying the thickness of the reactive sandpacks and pumping rates, three additional simulation cases were constructed. All other model parameters were not varied and are listed in Table 3. Fig. 6(a)–(d) show the simulation results of the four cases. Aqueous 90 Sr concentrations were plotted along the radial axis of the well, and the results are displayed at days 200, 400, 600, and 800. In general, after 800 days, the results suggest that the maximum sorption capacities will be reached or nearly reached for a sandpack thickness of 1 m (Fig. 6a and b), and that the concentrations in the well will be equivalent or nearly equivalent to the concentrations in the aquifer. Increasing the thickness of the sandpack to 2 m predicts substantially retarded 90 Sr breakthrough, especially for the case with a lower pumping rate of 1.5 × 10−2 m3 s−1 (Case 3). For Case 1, both kinetic models forecast early breakthrough, but also slightly longer attenuation of 90 Sr concentrations in comparison to the equilibrium case. These results indicate that for the same sorption capacity, kinetic limitations retain sorption capacity to a later stage of well operation. However, from a practical perspective these differences are not significant. The phenomenon is more pronounced for velocity-dependent kinetic sorption. In Case 2, due to the increased pumping rate, the maximum sorption capacity was nearly reached after 600 days. This result suggests that such a configuration may be inadequate to provide effective treatment for

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

a)Case1

693

b)Case2 2E-12

-1

Sr-90concentrations(molL )

Sr-90concentrations(molL -1 )

2E-12

d200(V-Dependent) d200(Kinetic) d200(Equilibrium) d400(V-Dependent) d400(Kinetic) d400(Equilibrium) d600(V-Dependent) d600(Kinetic) d600(Equilibrium) d800(V-Dependent) d800(Kinetic) d800(Equilibrium) Remediationtarget

1E-12

0

0

1

d200(V-Dependent) d200(Kinetic) d200(Equilibrium) d400(V-Dependent) d400(Kinetic) d400(Equilibrium) d600(V-Dependent) d600(Kinetic) d600(Equilibrium) Remediationtarget

1E-12

0

2

0

1

2

Distance(m)

Distance(m)

d)Case4

c)Case3

-1

-1

Sr-90concentrations(molL )

2E-12

Sr-90concentrations(molL )

2E-12

d200(V-Dependent) d200(Kinetic) d200(Equilibrium) d400(V-Dependent) d400(Kinetic) d400(Equilibrium) d600(V-Dependent) d600(Kinetic) d600(Equilibrium) d800(V-Dependent) d800(Kinetic) d800(Equilibrium) Remediationtarget

1E-12

0

0

1

2

3

4

Distance(m)

d200(V-Dependent) d200(Kinetic) d200(Equilibrium) d400(V-Dependent) d400(Kinetic) d400(Equilibrium) d600(V-Dependent) d600(Kinetic) d600(Equilibrium) d800(V-Dependent) d800(Kinetic) d800(Equilibrium) Remediationtarget

1E-12

0

0

1

2

3

4

Distance(m)

Fig. 6. Simulated aqueous 90 Sr concentration distributions along the distance from the pumping well screen using different pumping rates and reactive sandpack thicknesses.

extended time periods, independent of kinetic effects. In Case 3, the thickness of the reactive sandpacks was doubled to 2 m while maintaining the pumping rate. Therefore, after 800 days of simulation only partial breakthrough occurred regardless of the chosen model. In Case 4, both sandpack thickness and pumping rate were doubled and in this case the differences between the velocitydependent kinetic model and the other two models are largest. Due to the thickness of the sandpack, the kinetic rate (adopted from column C4) in the constant rate kinetic model provides simulation results more similar to the results of the equilibrium model. However, for the velocity-dependent kinetic model, aqueous 90 Sr is more rapidly transported toward the pumping well although clinoptilolite along the flow path still has the ability to retain the contaminant. Moreover, the kinetic effects increase, as velocity increases toward the well screen. This case clearly demonstrates that a thorough understanding of velocity-dependent sorption can be important for evaluating 90 Sr breakthrough with respect to a set target concentration. Fig. 7(a) and (b) show the simulated adsorbed 90 Sr concentrations for Case 3 where the sandpack had a radius of 2 m and the pumping rate was 1.5 × 10−2 m3 s−1 . The results indicate that near the well screen (Fig. 7(a)) the accumulation rates of the sorbed 90 Sr vary as a function of the sorption model; the velocity-dependent kinetic model predicts the highest accumulation rate. The major reason is that more 90 Sr mass was transported to the region near the well screen without being adsorbed within the outer layers

of the reactive sandpack. However, at a greater distance from the well center (Fig. 7(b)), sorbed 90 Sr masses obtained from the different models are almost identical, suggesting that kinetic limitations have a smaller effect on total adsorbed mass at lower pore-water flow rates. All four simulation cases indicate that installing reactive sandpacks will lead to a reduction of the 90 Sr contaminant levels for a period of time. It is important to determine how long the sandpacks will last to maintain 90 Sr concentrations in the extracted groundwater below a specified regulatory standard. For example, Health Canada suggests the maximum 90 Sr concentration in the safe drinking water to be 5 Bq L−1 (1.09 × 10−14 mol L−1 ). Here we define breakthrough when the 90 Sr in the pumping well exceeds that drinking water standard. Considering that the remediation target in a dewatering project might be set at a higher concentration, an alternative standard of 50 Bq L−1 was also used. Table 4 lists the simulated breakthrough times for the four cases and the differences between the velocity-dependent kinetic model and the other two models. In all cases, increasing the thickness of the sandpack and reducing the pumping rate increased the breakthrough times as expected. For example in the velocity-dependent model, if the pumping rate was 0.015 m3 s−1 , increasing the thickness of the sandpack from 1 m to 2 m increased the breakthrough times from 31 days to 502 days. If the regulatory standard were set at 50 Bq L−1 , the efficiency of the 90 Sr attenuation increased 300% in Case 1, 50% in Case 2, 59% in Case 3 and 1500% in Case 4. One important aspect

694

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695

a) 0.17 m (Case 3)

b) 1.13 m (Case 3)

3E-09

Sr-90 concentrations (mol L-1 )

Sr-90 concentrations (mol L-1 )

V-Dependent Kinetic Equilibrium

2E-09

1E-09

0

0

200

400

600

90

4E-08

2E-08

0

800

Days Fig. 7. Simulated sorbed 1.13 m.

6E-08

V-Dependent Kinetic Equilibrium 0

200

400

600

800

Days

Sr concentrations using different model formulations in Case 3 (sandpack radius: 2 m; pumping rate: 1.5 × 10−2 m3 s−1 ) at (a) 0.17 m and (b)

Table 4 Breakthrough time in days for the four field scale simulation cases.

Reactive sand pack thickness (m) Pumping rate (m3 s−1 ) Equilibrium model Kinetic model Velocity-dependent kinetic model a

5 Bq L−1 a 50 Bq L−1 5 Bq L−1 a 50 Bq L−1 5 Bq L−1 a 50 Bq L−1

Case 1

Case 2

Case 3

Case 4

1 1.5 × 10−2 142 215 106 192 31 124

1 3.0 × 10−2 71 108 40 86 0.04 0.06

2 1.5 × 10−2 713 >800 660 >800 502 >800

3 3.0 × 10−2 356 490 307 459 11 180

Breakthrough is taken at the time when the drinking water standard was reached 5 Bq L−1 (1.09 × 10−14 mol L−1 ).

of using the velocity-dependent model is that the breakthrough time decreased significantly to the order of hours in both regulatory standards when the pumping rate was 0.03 m3 s−1 and the sandpack thickness was 1 m (Case 2). Near instantaneous breakthrough in Case 2 is mainly due to two factors: (1) large amount of the 90 Sr mass was transported to the well screen due to the nonequilibrium sorption along the flow path; (2) the highest velocity near the well screen significantly decreased the kinetic rate and the retardation was minimized. Results suggest that the Case 3 design will be able to maintain groundwater quality at drinking water level for more than a year and below 50 Bq L−1 for over 2 years. In a dewatering project, the selection of the pumping rate and the sandpack thickness depend on a number of other factors such as local plume concentration, construction cost and duration of the dewatering project. If the project construction time is on the order of months, and if 50 Bq L−1 were set as the remediation target, the simulations suggest that three out of the four cases considered will be able to maintain groundwater quality at acceptable levels. Here the estimation of the breakthrough time is slightly conservative because radioactive decay will reduce the groundwater 90 Sr concentrations by 2% in one year and 7% in three years. Unless the time span of the dewatering is on the order of decades, which is unlikely, the impacts of decay can be neglected. Some uncertainties in the model calibration may affect the use of the model for field scale design simulations. For examples, the model is only valid for conditions at the Chalk River site because the calibrated Kd values are apparent partitioning coefficients. Variations in other chemical parameters, especially, in this instance, Ca2+ , certainly affect mobility of 90 Sr. If the local water chemistry changes significantly, distribution coefficients must be updated accordingly. Nevertheless, kinetic limitations and the power–law relationship

between the kinetic rates and pore water velocities agreed with previous observations. The present results suggest that the dependency of the kinetic rates on flow velocities must be considered in the design of reactive sandpacks. Previous studies [31–33] indicate that parameter upscaling is necessary for modeling contaminant transport at larger scales because the complexity of physical properties of the media (e.g., effective diffusion coefficients, retardation factor, porosity, and reactive mineral facies) increases. In general, if a statistical model of the media can be established, efficient parameters representing the “ensemble” conditions can be obtained through analytical solutions. In this study, considering that the clinoptilolite is relatively homogeneous, and the scale “leap” from the column (10−1 m) to the sandpack (101 m) is relatively small, upscaling was not considered. However, uncertainties associated with scaling effects should be taken into consideration when designing field applications of this technology. 5. Conclusions Granular clinoptilolite has a high selectivity toward Sr2+ and has proven efficient in attenuating a 90 Sr plume at Chalk River. Most of the previous studies using clinoptilolite were conducted at the batch scale. In this paper, we presented the results of four in situ column experiments conducted in an actual 90 Sr plume at the Chalk River Laboratories. Breakthrough curves were obtained for all four columns and non-competitive equilibrium and kinetic sorption models with apparent Kd were calibrated. The calibrated Kd values in the equilibrium model showed that Kd was affected by chemical conditions in the ambient groundwater (e.g., Ca2+ concentrations). The kinetic model captured the early breakthrough as well as the

J. Yin et al. / Journal of Hazardous Materials 280 (2014) 685–695 90 Sr

oscillations in the effluent. Moreover, calibrated kinetic rates were inversely correlated to the averaged pore water velocities. The short contact time between aqueous 90 Sr and clinoptilolite due to the higher pore water velocity and bypass of intra-granular exchange sites are the most likely reasons for the slower kinetic rates at high velocities. An experimentally derived regression equation, a velocity-dependent kinetic model was developed. All three sorption models were used to simulate 90 Sr attenuation for hypothetical dewatering well scenarios with the goal of evaluating the conditions under which reactive sandpacks can provide a suitable water treatment technology. The results suggest that kinetic effects, in particular near the well screen, results in early breakthrough. Therefore, field application of this technique may require thick sandpacks or slower pumping rates if near-complete contaminant removal for many months is required. Simulation results suggest that 90 Sr concentrations may approach the drinking water standard in one month for the base case with 1-m thick sandpacks at a pumping rate of 1.5 × 10−2 m3 s−1 . But for a remediation target of 50 Bq L−1 , the breakthrough time extended to 4 months. Increasing the sandpack thickness efficiently increases the breakthrough time, but would also substantially increase construction costs. Therefore, it is suggested that in an actual dewatering project, pre-modeling is an efficient way to provide quantitative and economic design parameters. Acknowledgements This research was supported by the Atomic Energy of Canada Limited, Chalk River Laboratories through research agreement F1103864. We also would like to thank Dr. Alan Rabideau for providing the batch experiment data and Karen Sharp for editorial improvements. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jhazmat. 2014.07.073. References [1] R.E. Jackson, K.J. Inch, The in-situ sorption of 90 Sr in a sand aquifer at the Chalk River Nuclear Laboratories, J. Contam. Hydrol. 4 (1989) 27–50. [2] R.W.D. Killey, J.H. Munch, Radiostrontium migration from a 1953-54 liquid release to a sand aquifer, Water Poll. Res. J. Canada 22 (1987) 107–128. [3] R.E. Jackson, K.J. Inch, R.J. Patterson, R.W.D. Killey, Field measurement and geochemical interpretation of 90 Sr sorption in a contaminated aquifer, EOS 65 (1984) 207. [4] R.E. Jackson, K.J. Inch, Partitioning of strontium-90 among aqueous and mineral species on a contaminated aquifer, Environ. Sci. Technol. 17 (1983) 231–237. [5] A.M. El-Kamash, M.R. El-Naggar, M.I. El-Dessouky, Immobilization of cesium and strontium radionuclides in zeolite-cement blends, J. Hazard. Mater. B136 (2006) 310–316. [6] R. Sureda, X. Martínez Lladó, M. Roviraa, J. de Pabloa, I. Casas, J. Giménez, Sorption of strontium on uranyl peroxide: implications for a high-level nuclear water repository, J. Hazard. Mater. 181 (2010) 881–885. [7] D.R. Lee, D.S. Hartwig, Interception of a groundwater plume containing strontium-90, in: Waste Management, Decommissioning and Environmental Restoration of Canada’s Nuclear Activities, Toronto, Ontario, Canada, September 11–14, 2011, Canadian Nuclear Society, 2011. [8] D.R. Lee, D.J.A. Smyth, S.G. Shikaze, R. Jowett, D.S. Hartwig, C. Milloy, Walland-curtain for passive collection/treatment of contaminant plumes, in: Designing and Applying Treatment Technologies: Remediation of Chlorinated and Recalcitrant Compounds. Proceedings of the First International Conference, Monterey, California, May 18–21, Battelle Press, Columbus, OH, 1998, pp. 77–84.

695

[9] J.P. Powers, A.B. Corwin, P.C. Schmall, W.E. Kaeck, Construction Dewatering and Groundwater Control: New Methods and Applications, 3rd ed., John Wiley & Sons, Hoboken, New Jersey, 2007, 638 pp. [10] E.H. Priebe, The Design of Reactive Sandpacks for the Attenuation of Radiostrontium in Groundwater for Dewatering Application (MSc Research Paper), Department of Earth and Environmental Sciences, University of Waterloo, Ontario, 2010. [11] D.R. Lee, D.S. Hartwig, Zeolite prevents discharge of strontium-90contaminated groundwater, in: Waste Management, Decommissioning and Environmental Restoration for Canada’s Nuclear Activities, Current Practices and Future Needs, Canadian Nuclear Society, Ottawa, Ontario, Canada, May 8–11, 2005. [12] E.H. Priebe, D.R. Lee, Reactive Sandpacks for In-situ Treatment of Construction Dewatering Effluent, AECL Nuclear Review 2 (2) (2013) 85–87. [13] N. Toride, F.J. Leij, M.Th. van Genuchten, The CXTFIT Code for Estimating Transport Parameters from Laboratory or Field Tracer Experiments: Version 2.0. U.S. Department of Agriculture. Research Report No. 137, 1995, pp. 121. [14] K.U. Mayer, E.O. Frind, D.W. Blowes, Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions, Water Resour. Res. 38 (2002), http://dx.doi.org/10.1029/2001WR000862. [15] C.L. Steefel, A.C. Lasaga, A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems, Am. J. Sci. 294 (1994) 529–592. [16] S.-W. Jeen, K.U. Mayer, R.W. Gillham, D.W. Blowes, Reactive transport modeling of trichloroethene treatment with declining reactivity of iron, Environ. Sci. Technol. 41 (2007) 1432–1438. [17] K.U. Mayer, S.G. Benner, D.W. Blowes, Process-based reactive transport modeling of a permeable reactive barrier for the treatment of mine drainage, J. Contam. Hydrol. 85 (2006) 195–211. [18] A.M. El-Kamash, Evaluation of zeolite A for the sorptive removal of Cs+ and Sr2+ ions from aqueous solutions using batch and fixed bed column operations, J. Hazard. Mater. 151 (2008) 432–445. [19] J. Doherty, PEST: Model-Independent Parameter Estimation, Watermark Numerical Computing, Brisbane, Australia, 2008. [20] A.J. Rabideau, J. Van Benschoten, A. Patel, K. Bandilla, Performance assessment of a zeolite treatment wall for removing 90 Sr from groundwater, J. Contam. Hydrol. 79 (2005) 1–24. [21] M. Fuhrmann, D. Aloysius, H. Zhou, Permeable, subsurface sorbent barrier for 90 Sr: laboratory studies of natural and synthetic materials, Waste Manage. 15 (1995) 485–493. [22] J. Samper, Z. Dai, J. Molinero, M. Garcia-Gutierrez, T. Missana, M. Mingarro, Inverse modeling of tracer experiments in FEBEX compacted Ca-bentonite, Phys. Chem. Earth 31 (10) (2006) 640–648. [23] W. Dong, S.C. Brooks, Determination of the formation constants of ternary complexes of uranyl and carbonate with alkaline earth metal (Mg2+ , Ca2+ , Sr2+ , and Ba2+ ) using anion exchange method, Environ. Sci. Technol. 40 (2006) 4689–4695. [24] J. Yin, R. Haggerty, D.L. Stoliker, D.B. Kent, J.D. Istok, J. Greskowiak, J.M. Zachara, Transient groundwater chemistry near a river: effects on U(VI) transport in laboratory column experiments, Water Resour. Res. 47 (2011) W04502, http://dx.doi.org/10.1029/2010WR009369. [25] Z. Dai, A. Wolfsberg, P. Reimus, H. Deng, E. Kwicklis, M. Ding, D. Ware, M. Ye, Identification of sorption processes and parameters for radionuclide transport in fractured rock, J. Hydrol. 414–415 (2012) 220–230. [26] X. Zhang, J.W. Crawford, I.M. Young, Does pore water velocity affect the reaction rates of adsorptive solute transport in soils? Demonstration with pore-scale modeling, Adv. Water Res. 3 (2008) 425–437. [27] S. Kim, H. Ha, N. Choi, D. Kim, Influence of flow rate and organic carbon content on benzene transport in a sandy soil, Hydrol. Process. 20 (2006) 4307–4316. [28] M.A. Maraqa, Prediction of mass-transfer coefficient for solute transport in porous media, J. Contam. Hydrol. 50 (2001) 1–19. [29] Lappala E.G., Healy R.W., Weeks E.P., Documentation of Computer Program VS2D to Solve the Equations of Fluid Flow in Variable Saturated Porous Media. Water Resources Investigation Report 83-4099, U.S. Geological Survey, Denver, CO, 1987. [30] T.H. Henderson, Numerical Modeling of Density-Driven Chemical Oxidation of Chlorinated Solvents (Ph.D. Thesis), University of British Columbia, Vancouver, Canada, 2009. [31] Z. Dai, A. Wolfsberg, Z. Lu, H. Deng, Scale dependence of sorption coefficients for contaminant transport in saturated fractured rock, Geophys. Res. Lett. 36 (2009), http://dx.doi.org/10.1029/2008GL036516. [32] H. Deng, Z. Dai, A. Wolfsberg, Z. Lu, M. Ye, P. Reimus, Upscaling of reactive mass transport in fractured rocks with multimodal reactive mineral facies, Water Resour. Res. 46 (2010) W06501, http://dx.doi.org/10.1029/2009WR008363. [33] H. Deng, Z. Dai, A.V. Wolfsberg, M. Ye, P. Stauffer, Z. Lu, E. Kwicklis, Upscaling retardation factor in hierarchical porous media with multimodal reactive mineral facies, Chemosphere 91 (2013) 248–257.