Groundwater transport modeling with nonlinear sorption and intraparticle diffusion

Accepted Manuscript Groundwater transport modeling with nonlinear sorption and intraparticle diffusion Anshuman Singh, Richelle M. Allen-King, Alan J. Rabideau PII: DOI: Reference:

S0309-1708(14)00077-3 http://dx.doi.org/10.1016/j.advwatres.2014.04.010 ADWR 2191

To appear in:

Advances in Water Resources

Received Date: Revised Date: Accepted Date:

30 September 2013 6 April 2014 8 April 2014

Please cite this article as: Singh, A., Allen-King, R.M., Rabideau, A.J., Groundwater transport modeling with nonlinear sorption and intraparticle diffusion, Advances in Water Resources (2014), doi: http://dx.doi.org/10.1016/ j.advwatres.2014.04.010

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Groundwater transport modeling with nonlinear sorption and intraparticle diffusion Anshuman Singha,, Richelle M. Allen-Kingb, Alan J. Rabideaua*

a

Department of Civil, Structural and Environmental Engineering, University at Buffalo, Buffalo,

NY 14260.

b

Department of Geology, University at Buffalo, Buffalo, NY 14260.

*Corresponding Author

email addresses: [email protected]; [email protected]; [email protected]

ABSTRACT Despite recent advances in the mechanistic understanding of sorption in groundwater systems, most contaminant transport models provide limited support for nonideal sorption processes such as nonlinear isotherms and/or diffusion-limited sorption. However, recent developments in the conceptualization of “dual mode” sorption for hydrophobic organic contaminants have provided more realistic and mechanistically sound alternatives to the commonly used Langmuir and Freundlich models. To support the inclusion of both nonlinear and diffusion-limited sorption processes in groundwater transport models, this paper presents two numerical algorithms based on the split operator approach. For the nonlinear equilibrium scenario, the commonly used twostep split operator algorithm has been modified to provide a more robust treatment of complex 1

multi-parameter isotherms such as the Polanyi-partitioning model. For diffusion-limited sorption, a flexible three step split-operator procedure is presented to simulate intraparticle diffusion in multiple spherical particles with different sizes and nonlinear isotherms. Numerical experiments confirmed the accuracy of both algorithms for several candidate isotherms. However, the primary advantages of the algorithms are: (1) flexibility to accommodate any isotherm equation including “dual mode” and similar expressions, and (2) ease of adapting existing grid-based transport models of any dimensionality to include nonlinear sorption and/or intraparticle diffusion. Comparisons are developed for one-dimensional transport scenarios with different isotherms and particle configurations. Illustrative results highlight (1) the potential influence of isotherm model selection on solute transport predictions, and (2) the combined effects of intraparticle diffusion and nonlinear sorption on the plume transport and flushing for both singleparticle and multi-particle scenarios. 1 INTRODUCTION Of the processes that affect the migration of dissolved contaminants in the subsurface, sorption to aquifer solids is believed to play a key role, especially in the absence of nonaqueous phase contamination. In particular, the influence of sorption on aquifer restoration has long been recognized (e.g. [1-5]), and the release of sorbed contaminants is an important component of modern multi-compartment conceptualizations of remediation processes (e.g. [6]). The prediction of cleanup times for hydrophobic organic contaminants (HOCs) can be strongly affected by assumptions used to model the sorption process, which typically include specification of an equilibrium sorption isotherm and, if applicable, a rate model [7-10]. Although potentially influential aspects of sorption phenomena such as nonlinear equilibrium and/or kinetic limitations have been conceptualized for subsurface systems [11-15], only a few

2

published field scale modeling studies have incorporated these processes [9, 16, 17], although some public domain software packages support a subset of the available models [18]. As the understanding of HOC sorption has evolved, recent developments have introduced complex isotherm expressions involving multi-parameter “dual mode” or “distributed reactivity” mechanisms [19-21]. However, in most groundwater contaminant transport models, HOC sorption is represented as a linear process, or by a limited number of nonlinear isotherms such as the Freundlich [22] or Langmuir [23] equation [18, 24-27]. Also, the combination of intraparticle diffusion and nonlinear isotherms, while recognized as an appealing mechanistic conceptualization of “nonideal” sorption processes (e.g.[13, 28-31]) has been modeled primarily at the laboratory scale [31-35], with only a few studies conducted for intraparticle diffusion in two-dimensional domains [9].

The disconnect between the current mechanistic understanding of subsurface sorption and the lack of relevant modeling studies is most likely attributable to the significant computational burden of coupling field scale solute transport with grain scale diffusion processes. Thus, the primary purpose of this work is to present a modified computational algorithm designed to incorporate both aspects of “nonideal” sorption in a manner that can be scaled in a straightforward fashion to two- and three-dimensional applications. Although the focus of this work is on demonstrating the concepts through one-dimensional scenarios, the numerical algorithms have been designed to facilitate eventual extension to multidimensional field scale transport, for which the relative importance of grain scale processes remains largely unexplored.

To represent nonequilibrium (slow) sorption, we have adopted a traditional pore diffusion formulation, as detailed in Methods, which assumes that the rate of the sorption process is 3

controlled by the physical transport of aqueous contaminants to sorbents located within the interior of soil grains, where local equilibrium is achieved. However, for nonlinear equilibrium, we have extended standard transport formulations to incorporate the emerging class of dual mode isotherms. As reviewed by Allen-King et al. [21], these models combine a nonlinear term to represent surface adsorption to “hard” carbon with a linear term based on partitioning to “soft” carbon. In particular, the Polanyi-partitioning (PP) isotherm presented below in Eq. (1) has been emphasized because of its mechanistically grounded “pore filling” expression for surface adsorption component [20, 36].

qs = Q p 10

    

     



S   A log w   Cl  

B (1) + K P Cl

where qs is the mass fraction of the contaminant in the sorbent phase [M/M], Cl is the aqueous concentration of solute [M/L3], Qp and A, B are Polanyi parameters, Sw is the solute aqueous solubility [M/L3], and KP is the partition coefficient used in the PP formulation [L3/M] (see [21] and [36]).

At the time of this writing, the PP equation has not been utilized in published groundwater transport models and is therefore the focus of the numerical developments associated with this work. For comparison purposes, we also consider the commonly used nonlinear Freundlich and Langmuir isotherm expressions:

Freundlich:

qs = K f Clnf

(2)

4

Langmuir:

qs =

Q 0C l 1 + bCl

(3)

where Kf and nf are the Freundlich parameters and Q0 and b are the Langmuir parameters [L3/M]. It is noted that both isotherms can be reduced to the more commonly used linear isotherm (qs = Kd Cl) with appropriate specification of parameters.

In the next section, we present an efficient and accurate split-operator algorithm that can incorporate any isotherm expression into a intraparticle diffusion formulation, with a modular structure that can be coupled in a straightforward fashion to existing transport codes. The threestep numerical algorithm is based on the common split-operator decoupling of macroscopic advective-dispersive transport from local reaction processes, extended via a similar approach to decouple the local diffusion process from nonlinear sorption at intraparticle computational nodes. This approach avoids the use of ODE solvers and accommodates nonlinear isotherms by solving a single nonlinear algebraic mass balance equation for each macroscopic spatial node (for equilibrium sorption) or at intraparticle nodes (for diffusion limited sorption). Advantages of this approach are similar to other split-operator methods for reactive transport modeling, including reduced demand for computational processing time and memory, suitability for parallel processing, ease of code maintainability, and flexibility to incorporate new combinations of reactions [37]. The three-step algorithm is demonstrated for one-dimensional solute transport problems incorporating the nonlinear Freundlich and the dual mode PP isotherms, as implemented in public domain software [38].

5

2 MODEL FORMULATION 2.1 Governing Equations

The one-dimensional form of the advective-dispersive-reactive equation (ADRE), with sorption considered as the dominant “reaction,” can be represented as:

∂Cl ∂C ∂ 2C 1 − n  ∂Cs  = −v l + D 2 l −   ∂t ∂x ∂x n  ∂t 

(4)

where Cl and Cs are the solute concentrations in the aqueous and sorbed phases, respectively [M/L3], each averaged with respect to the respective phase volume (aqueous or solid); v is the groundwater velocity [L/T]; D is the dispersion coefficient [L2/T]; t is time [T]; x is distance from the domain entrance [L]; and n is the aquifer porosity. The formulation can readily be extended for transport in two- or three-dimensions without modifying the sorption term. At equilibrium, the mass fraction of sorbed contaminant (qs) is a function of aqueous concentration (Cl) as expressed by the governing isotherm (equations (1)-(3)). The volume averaged sorbed phase concentration (Cs) in Eq. (4) is related to qs through the relationship ρb qs = (1 - n) Cs, where ρb is the bulk density of the aquifer material [M/L3].

When sorption is rate limited, the process can be represented as aqueous diffusion into intraparticle pore fluid coupled with rapid sorption to organic matter (both hard and soft, if applicable) located on the adjacent intraparticle pore walls. For an assemblage of spherical particles of different sizes and sorption characteristics, the governing equations can be expressed as follows [13, 39-42]:

6

 ∂2c ∂clk 2 ∂clk = Dak  2lk + ∂t rk ∂rk  ∂rk

 ρ bk ∂qsk  −  n pk ∂t

(5)

where the subscript k represents a class of particles with specified diffusion and sorption characteristics; rk denotes the radial distance from the particle center [L]; ρbk is the particle bulk density [M/L3]; clk is the solute concentration in the intraparticle fluid [M/L3]; qsk is the local sorbed mass fraction [M/M]; Dak is the effective solute intraparticle diffusion coefficient (including tortuosity effects) [L2/T] and npk is the intraparticle porosity. As before, the equilibrium relationship between the local aqueous concentration clk and sorbed fraction qsk is expressed using one of the isotherms given in Eqs. (1) – (3). The volume averaged macroscopic sorbed concentration Cs is related to the intraparticle concentrations as follows:

  ρ χ b k Cs = ∑  4 3  k  π Rk ρ bk 3

∫

Rk

0

  (n pk clk + ρ bk q sk )4πrk2 drk    

(6)

where, Rk is the particle radius for class k [L]; χk is the mass fraction of the sorbent, and ρb is the macroscopic bulk density [M/L]. The number of particles per unit volume, Nk, is given by:

Nk =

ρb χ i 4 3 πRk ρ bk 3

2.2 Solution Methods

For equilibrium (fast) sorption, the incorporation of nonlinear isotherm expressions into the advective-dispersive-reactive-equation (ADRE) results in a system of nonlinear equations that

7

(7)

must be solved numerically. The most common computational approach is to modify the advection/dispersion terms using a concentration-dependent retardation factor. Although iteration or frequent updating of the retardation factor is needed for accurate simulation of highly nonlinear systems, the approach has been successfully demonstrated for the Langmuir and Freundlich isotherms [e.g., 18, 43]. An alternative numerical strategy for modeling nonlinear sorption is to split the transport and reaction (sorption) operators, which generates smaller subsystems of nonlinear nodal sorption equations that are particularly suitable for parallel computing environments [44]. For rapid sorption, the process can be modeled as rate-limited, using a large sorption rate constant to mimic equilibrium conditions. Such an approach can introduce an additional programming burden, but the solution of the sorption subsystem is readily accomplished using standard ordinary differential equation (ODE) solvers e.g. [24 – 26]. Such split-operator schemes for the solution of the ADRE have been well established and studied for a variety of reaction scenarios [44 – 48]. However, the use of such “fast kinetic” approximations to simulate nonlinear sorption has two drawbacks: (1) the overhead associated with the use of an ODE solver for a relatively simple system of two nodal equations, and (2) the need to specify a fictitious rate constant that is large enough to ensure equilibrium behavior but not so large as to require excessive numerical iteration in the ODE solution.

The split-operator approach described above can be used to model rate limited sorption if the sorbed phase is considered as a single compartment [e.g., 18]. However, a more complex mathematical formulation is required if the rate limiting process is treated as a diffusion process (Eq. 5), resulting in considerable additional computation. Consequently, most previous work on coupling contaminant transport with intraparticle diffusion has emphasized one-dimensional 8

macroscale domains and simplified isotherms, as reviewed in [44]. Exceptions include the work of Rabideau and Miller [9] who studied two-dimensional radial transport, including some applications with nonlinear isotherms, and the work of Liedl and Ptak [40], who considered onedimensional transport for three particles of equal sizes and different sorption characteristics represented by Freundlich isotherms. Both applications utilized a split-operator numerical approach that could be readily extended to 3D transport. However, a realistic field scale simulation that includes intraparticle diffusion in multiple particle classes of different physical and sorptive properties would typically require excessive computational resources using existing methods and has not been reported.

2.2 Split-operator formulations

In this work, we extend previously developed split-operator methods to provide greater flexibility in accommodating both a diffusion-based kinetic formulation and nonlinear sorption isotherms. Eq. (4) can be split into two sequential steps, expressed by operators L1(Cl) for advection-dispersion and L2(Cl) for sorption, as given by:

∂Cl ∂C l ∂ 2 Cl = L1 (Cl ) = −v +D 2 ∂t ∂x ∂x

(8)

∂Cl 1 − n ∂Cs = L2 (Cl ) = − n ∂t ∂t

(9)

In the two-step split-operator procedure, the advective-dispersive portion of ADRE in Eq. (8) is solved over one time step ∆t using any suitable numerical procedure, yielding an intermediate

9

concentration vector C*. For this work, a Galerkin finite element numerical scheme was used with uniform spatial discretization of ∆x, combined with a Crank-Nicholson finite difference approximation of the time derivative. For the next sequential step (L2), the reaction equation for nonideal sorption is solved separately at each spatial node over one timestep (∆t). The formulation of Eq. (9) for equilibrium (fast) sorption depends on the isotherm relationship given in Eqs. (1) – (3). For this class of problems, Eq. (9) can be replaced by a mass balance at each spatial node i:

nC lT,1i +1 + ρ b q sT1,i+1 = nC l*,i + ρ b q sT,1i

(10)

where the superscript T1 and T1 + 1 denote the beginning and end of the ADRE time step, respectively; ClT,i +1 and qsT,i+1 are the aqueous phase concentration and sorbed mass fraction at time 1

1

T1+1; q sT,1i is sorbed mass fraction at time T1; and Ci* is the intermediate aqueous concentration

obtained by solving the transport operator L1. Using the above mass balance, Eq. (9) reduces to a single nonlinear nodal equation depending on the form of the isotherm expression, which is used to replace qsT,i+1 with the appropriate function of ClT,i +1 . 1

1

For nonequilbrium sorption modeled as intraparticle diffusion, Eq. (9) is represented by the radial diffusion equations, which are further split into two sequential steps expressed by operators L3(clk) for diffusion and L4(clk) for nonlinear sorption:

 ∂ 2 clk 2 ∂clk  ∂clk  = L3 (clk ) = Dak  2 + ∂t r ∂ r ∂ r k k   k

10

(11)

∂clk ρ ∂q = L4 (clk ) = − bk sk ∂t n pk ∂t

(12)

For this work, a finite difference numerical solution to Eq. (11) was used (Appendix A.2), which requires discretization of each particle into radial nodes. Following a similar strategy to the nonlinear equilibrium case, Eq. (12) can be represented by a mass balance at each internal node j of the spherical particle:

n pk clk* , j + ρ bk q skT2 , j = n pk clkT2,+j1 + ρ bk q skT2 ,+j1

(13)

where the superscript T2 and T2+1 represent the times at the beginning and end of the discretized diffusion time sequence (obtained by subdividing one ADRE time interval); clkT ,+j1 and qskT ,+j1 are 2

2

the aqueous phase concentration and sorbed mass fraction from diffusion-reaction time T2+1 at particle node j; q skT , j is the sorbed mass fraction from diffusion-reaction time T2; and clk* , j is 2

obtained by solving the operator L3. Similar to the previous case, Eq. (12) reduces to a nonlinear algebraic equation depending on the isotherm expression in Eqs. (1) – (3). Analysis of test problems is aided by defining the following dimensionless groupings: Pe =

v∆x D

(14a)

Cr =

v ∆t ∆x

(14b)

Da =

Da L vR 2

(14c)

D a ∆t d ∆r 2

(14c)

F=

11

Where Pe is the dimensionless mesh Peclet number, ∆x is the grid spacing in the macroscopic (advective-dispersive) domain [L], Cr is the dimensionless Courant number, ∆t is the transport time step [L], Da is the dimensionless Damkohler number for intraparticle diffusion, L is the advection length scale (one-dimensional domain length for the examples of this work), R is the particle radius [L], F is the dimensionless Fourier number, ∆td is the intraparticle diffusion time step (subdivided transport time step) [T], and ∆r is the intraparticle node spacing [L].

In general, the standard Galerkin FEM formulation used for the transport solution, as well as the FD formulations used in comparison models, require Pe < 2 for numerical accuracy. In addition, the explicit finite difference formulation used to solve the intraparticle diffusion equation (Appendix A.2) imposes a stability requirement on the subdivided time step of F < 0.5. Although the FEM and FD algorithms require Cr < 1 to maintain accuracy, the control of operator-splitting errors typically necessitate a smaller time step, and the influence of Cr on problems of interest to this work is explored in the subsequent examples. The Damkohler number is typically used to characterize the degree of nonequilibrium caused by intraparticle diffusion, with D a ≈ 100 or larger commonly associated with near-equilibrium conditions [e.g., 13].

3 SOLUTION METHODS AND TESTING 3.1 Equilibrium sorption

Each of the split operator algorithms has been implemented in the public domain MOUSER software, which uses a flexible quadratic FEM solution for the advective-dispersive transport operator [38]. For nonlinear equilibrium sorption, Eq. (8) is first solved to yield the intermediate Cl* at each node of the discretized domain. The nonlinear nodal mass balance 12

equations (10) are then solved independently for ClT,i +1 as a function of the known values Cl* and 1

q sT,1i (with q lT,2i +1 replaced by the isotherm equation), using the bisection method outlined in

Appendix A.1. For this work, a constant transport time step (∆t) was used for both equilibrium and intraparticle diffusion simulations. However, the numerical accuracy and computational efficiency of the split-operator algorithm could be further improved through the use of an adaptive time-stepping procedure as described by Gasda et al. [49].

Because closed-form solutions for the ADRE are unavailable for scenarios involving nonlinear isotherms, a variety of procedures were used to evaluate the accuracy and efficiency of the twostep algorithm, which is subsequently referred to as the mass balance (MB) procedure. First, nonlinear isotherm parameters were configured to mimic equivalent linear cases, and solutions were compared with analytical solutions for transport with linear sorption and idealized boundary conditions [50]. Second, for nonlinear cases, high resolution solutions with very small spatial and temporal discretization were developed to verify convergence of the two-step solution with mesh refinement. Third, internal mass balance calculations were implemented in the software to verify numerical mass conservation. Fourth, the MB algorithm was configured to replicate the Freundlich and Langmuir test problems developed by Grove and Stollenwork [51], as described in the MT3D software manual [18].

To illustrate the algorithm for a more realistic problem of interest, a one-dimensional test problem was developed for a step source, using isotherm parameters for aquifer material from Gull River [52], which were fit using the regression algorithms implemented in the ISOFIT program [53] with the data weighted proportional to errors estimated via eror propagation. 13

Results for all three isotherms were generated using the MB algorithm of the MOUSER program, with comparative results for the commonly used Freundlich and Langmuir isotherms provided by the public domain codes MT3D [18] and NIGHTHAWK [25]. For MT3D, an equivalent three-dimensional domain was configured to match the test problem, and both the total variation diminishing (TVD) and implicit finite difference (FD) methods were used to solve the ADRE. Both MT3D approaches resolve the isotherm nonlinearity using a lagged retardation factor that is updated during the iterative solution of the linearized ADRE.

In contrast,

NIGHTHAWK solves the one-dimensional ADRE using a sequential (non-iterative) splitoperator approach using the finite difference method for the transport operator and a “fast kinetic” approximation for the reactions. All comparisons utilized consistent spatial and temporal discretization.

3.2 Intraparticle diffusion with single particle type

When sorption is rate limited by intraparticle diffusion, the second step of the split-operator scheme requires solution of the nodal diffusion equations (11). In this work, an explicit finite difference scheme was developed, as summarized in Appendix A.2. Although the potential for parallel processing motivated the selection of an explicit algorithm to simulate intraparticle diffusion, the most effective parallelization strategy was to assign the entire reaction solution for each macroscopic node (or set of nodes) to a different processor. Because the solution of the diffusion equations accounted for the vast majority of computational time (typically > 90%), this simple scheme yielded near-linear speedup on a distributed memory workstation cluster for the test problems described in this work.

14

Spatial discretization of the nodal particles requires the solution of sets of coupled differential equations for each macroscopic node and transport time step (∆t), which is further subdivided into sorption time steps (∆td). The initial concentration in the bulk solution external to the particle is set as C* obtained from the solution of Eq. (8), while the initial concentrations and sorbed fractions at the discretized locations inside particles are obtained from the particle status at the end of the previous transport step. If the sorption isotherm is nonlinear, a third step is required to adjust intraparticle concentrations for sorption using Eq. (13), applying the bisection method as described previously. The method can incorporate any sorption isotherm into Eq. (13), and is demonstrated here for the Freundlich and PP isotherms.

The one-dimensional test problem developed for equilibrium sorption (Table 1) was extended to include intraparticle diffusion with particles of three different sizes (Table 2). These comparisons were conducted using consistent spatial and temporal discretization for the transport solution, and the intraparticle discretization was set such that the local mesh Fourier number based on the smallest radial spacing was always less than 0.5.

3.3 Intraparticle diffusion with mixture of particle types

A test case for a multi-particle scenario was developed based on the problem described in Liedl and Ptak [40], with Freundlich isotherm parameters for phenanthrene determined for laboratory mixtures of three soil types: sandstone, light limestone and dark limestone. This published test problem assigned equal proportions of the three types of particle. However, in the updated test problem presented here, the mass fractions of the particles were modified such that the three types of particles contained approximately equal initial sorbed masses when equilibrated with an initial aqueous concentration of 40 µg/L, which was applied across the entire macroscopic 15

domain. The model was then operated in “flushing” mode with a zero-concentration entrance boundary condition to mimic conditions associated with aquifer restoration. The diffusion coefficient for the dark limestone was also adjusted to explore its effect on the overall rate of mass removal from the system. The domain was divided such that Pe and Cr were 1.1 and 0.42 respectively, and the mesh Fourier number for the intraparticle diffusion model was 0.48, which necessitated 146 intraparticle diffusion time steps for each transport time step. For all of the multiparticle scenarios, the possible relationship between particle size and macrodispersion was neglected for clarity; a consistent macroscopic dispersion coefficient (D) was used for all simulations.

16

4 RESULTS 4.1 Preliminary simulations

Although few published solutions are available for the problems of interest to this work, the accuracy of the nonlinear equilibrium mass balance (MB) algorithm was assessed by comparisons against MT3D and NIGHTHAWK solutions for the one dimensional test problem of Grove and Stollenwerk [51] for Freundlich and Langmuir isotherms. The domain, source condition, and spatial/temporal discretizations were configured based on the example described in the MT3D manual [18]. For these problems, the spatial concentration distributions for a pulse source predicted by the MB algorithm were in good agreement with NIGHTHAWK and the MT3D (TVD) solution (Figure 1).

Figure 1: Predicted solute concentrations at the domain exit for the test problem of Grove and Stollenwerk [51] (Pe = 0.16 and Cr = 0.5).

17

Additional comparisons of the MB algorithm against MT3D and Nighthawk were performed for the test problem described in Table 1. For all of the test problems, the initial concentration in the domain was assumed to be zero and the boundary condition at the inlet and outlet were specified as Danckwerts (third type) and zero gradient, respectively. For the series of related test problems described in this and subsequent sections, split-operator convergence of the MB algorithm was established heuristically by repeated simulations using smaller time steps until simulated concentration distributions were graphically identical. Four examples are shown in Figure 2, including a intraparticle diffusion scenario using the additional parameters summarized in Table 2.

Although the effect of time step on accuracy was slightly different for each case, graphical agreement was observed between Cr = 0.05 and Cr = 0.1 for all problems discussed in this work. Similar relationships between operator-splitting errors and time step size have been reported for a variety of advection-dispersion-reaction problems (e.g. [44, 47, 54]). While a fairly consistent convergence threshold was observed for a variety of test problems based on the Table 1 sorption parameters, a focused heuristic analysis or adaptive error control [49] would be necessary to establish the appropriate temporal discretization for other problems, particularly scenarios involving greater sorption capacity and/or highly nonlinear isotherms.

18

Figure 2: Simulated temporal concentration profiles at the domain exit under conditions summarized in Tables 1-2 for different Cr values. Variations on the test problem were developed to provide comparisons with the MT3D and Nighthawk models, both of which support the Freundlich and Langmuir isotherms. For these comparisons, a consistent spatial discretization was maintained (Pe = 0.56), while Cr was varied from 0.05 to 1.0. As shown in Figure 3 (Freundlich) and Figure 4 (Langmuir), comparisons of 19

the simulated spatial concentration distribution after 600 days of transport indicated close agreement among all solutions at Cr = 0.1 or smaller, except for the MT3D finite difference (FD) solution, which exhibited discernible numerical dispersion. All of the numerical solutions exhibited increasing separation from the high-resolution solution with increasing Courant numbers (Cr > 0.1) for both Freundlich and Langmuir isotherms. Although formal benchmarking was not performed, the computational times required for MT3D and Mouser were comparable and typically lower when compared to NIGHTHAWK, which is attributed to the use of a generic ODE solver for the fast kinetic approach. The mass balance errors reported by Mouser and MT3D were consistently on the order of 0.1% (mass balance errors are not reported by NIGHTHAWK).

20

Figure 3: Simulated 600-day spatial concentration distributions for Freundlich isotherm. The black dotted line represents the “high resolution” MB output at Cr = 0.1

21

Figure 4: Simulated 600-day spatial concentration distributions for Langmuir isotherm. The black dotted line represents the “high resolution” MB output at Cr = 0.1.

22

4.3 Dual mode sorption

The MB algorithm was further tested for the conditions of Table 1 using the PP isotherm, which had previously not been implemented in published contaminant transport models. For the test problem, the numerical solution of Eq. 10 exhibited slow or non-convergence using NewtonRaphson iteration, and the bisection procedure (Appendix A.1) was therefore implemented to guarantee convergence. The expected relationship between time step and accuracy was observed, as shown in Figure 5, which also compares simulation results for the Freundlich and Langmuir isotherms.

Figure 5: Simulated 600-day spatial concentration distributions for PP isotherm at different Courant numbers using MB algorithm (black lines) for equilibrium sorption. The grey solid line represents the corresponding spatial plot for the Freundlich isotherm, while the grey dotted line represents the corresponding spatial plot for the Langmuir isotherm. .

23

The ability to accommodate a variety of isotherms expressions allows for a focused examination of the isotherm selection process. Figure 6 shows predicted concentration-versus-time plots using the three isotherms, each fit from the same data [52] using a similar weighted least-squares hybrid regression procedure [53]. The differences in model predictions for isotherms obtained from the same laboratory data reflect, in part, the low concentration of the pulse source (30 µg/L) relative to the concentration range of the laboratory isotherm study, which emphasized higher aqueous concentrations. However, similar simulations using a higher-concentration pulse source (not shown) also generated discernable, albeit less dramatic, differences in predicted concentration plots. Although the primary purpose of this example was to demonstrate the ability of the algorithm to accommodate the four-parameter PP isotherm expression, additional investigations are ongoing to examine potential difference in transport predictions using alternative isotherm models parameterized from the same data.

24

Figure 6: Predicted concentrations at the domain exit under equilibrium sorption for different isotherms (Fr=Freundlich; PP= Polanyi-partitioning; Lang=Langmuir) with Cr = 0.5.

4.4 Intraparticle diffusion: single particle

The effect of intraparticle diffusion was examined using a series of simulations implemented for the conditions summarized in Table 3. For single-particle scenarios, Figure 7 (PP) and Figure 8 (Freundlich) illustrate the effect of particle size on the predicted temporal concentration profiles at a distance of 2.275 meters from the domain entrance. For the conditions of these simulations, the kinetic effects of intraparticle diffusion for the smallest size particle (0.7 mm) were minimal, corresponding to near-equilibrium behavior for both isotherms ( Da = 56) . In Figure 9, the spatial concentration distributions after 600 days of transport are shown for both isotherms and several particle sizes. As with the temporal simulations, the effect of intraparticle diffusion was 25

very small for r = 0.7 mm but became more discernable with increasing particle size. In general, the relative importance of intraparticle diffusion depends on several factors, including the macroscopic velocity, isotherm characteristics, and intraparticle porosity and tortuosity, which can vary for different materials for the same grain size and sorption characteristics. The comparisons of Figure 9 also highlight the potential role of isotherm selection, as the PP and Freundlich parameters were fit to the same data using consistent methods.

Figure 7: Predicted concentration at domain exit for intraparticle intraparticle diffusion, three particle sizes, and PP isotherm (Cr = 0.5).

26

Figure 8: Predicted concentration at domain exit with intraparticle diffusion, three particle sizes, and Freundlich isotherm (Cr = 0.5).

27

Figure 9: Simulated 600-day spatial concentration profiles for Cr = 0.1. The thick solid grey lines represent equilibrium (EQ) sorption without intraparticle diffusion; the thin grey lines represent intraparticle diffusion and the Freundlich (Fr) isotherm. The thick solid black line represents equilibrium sorption without intraparticle diffusion and the PP isotherm; the thin black lines represent solutions with intraparticle diffusion. 4.4 Intraparticle diffusion: mixture of particle types

Figure 10 and Figure 11 illustrate a “flushing” scenario in which an initially contaminated region of aquifer is flooded with zero-concentration water, similar to conditions associated with aquifer remediation, for the conditions summarized in Table 4. For the case in which a very low effective diffusion coefficient was assigned to the “dark limestone,” most of the sorbed mass remained in the “dark” particle while the other two particles released mass at a much faster rate, as shown by the black lines. The slow release of contaminant was influenced by both the highly nonlinear Freundlich isotherm and the low diffusion coefficient associated with dark limestone. In a related numerical experiment shown by the grey lines, the effective diffusion coefficient of 28

the dark limestone was increased and equated to that of the light limestone (e.g., as reflective of a different tortuosity factor). This modification resulted in a higher overall mass release for the mixture, even though the properties of the light limestone and sandstone were unaltered.

The multiparticle results indicate how differences in individual particle properties can affect the overall desorption behavior of the system in a manner that is not straightforward to predict. In both cases, contaminant removal was dominated by a small fraction of the highly nonlinear/slowly desorbing particles. The availability of suitable modeling tools will support future research to better understand the factors that limit mass removal during pump-and-treat and similar technologies, including the presence of small quantities of “super sorbents”. Such insights, for example, could lead to modified site investigation and/or remediation strategies that emphasize a highly influential component of aquifer materials.

29

Figure 10: Predicted concentration profiles at x = 0.2 m under “flushing” conditions and desorption from mixture of particles with intraparticle diffusion. Dak values are for dark limestone.

30

Figure 11: Fraction of total sorbed mass remaining in domain under “flushing” conditions and desorption from a mixture of particles influenced by intraparticle diffusion. The black sets of lines indicate a lower diffusion coefficient for the dark particle. The solid lines represent the total sorbed fraction in domain (across all sorbents). The uniform dashed lines represent the sorbed mass for sandstone, dash- dot lines represent light limestone, and dash- dot-dot lines represent the dark limestone.

5 Conclusions, limitations, and extensions

Updated split-operator numerical algorithms were developed and evaluated to solve the ADRE with nonlinear sorption and/or intraparticle diffusion, including dual mode isotherms and/or multiple particle types having different physical properties and sorption characteristics. Implementation of the two-step and three-step numerical solutions with suitable attention to temporal discretization temporal discretization (typically Cr = 0.1 – 0.5), consistently gave solutions that were close to more refined solutions (Cr ~ 0.05), with the expected improvement in accuracy as the Courant number was decreased. For scenarios involving intraparticle diffusion, numerical accuracy was also dependent on the number of radial nodes; a higher 31

number of radial nodes were typically required for highly sorbing particles (on the order of 20100 nodes, with higher number of nodes for larger particles). For realistic test problems of interest, the kinetic influence of intraparticle diffusion was increasingly evident as the particle radius was increased beyond 0.7 mm, although the equilibrium threshold will depend upon the macroscopic fluid velocity and diffusion/sorption characteristics.

Of potential relevance to aquifer restoration, the model demonstrated the effects of diffusion and sorption for mixtures of particles with different physical and sorption properties. The examples highlight the potential influence of a small percentage of highly sorbing particles on overall mass removal from the porous media system. Thus, future work based on the intraparticle diffusion model could potentially contribute to an improved understanding of the contaminant behavior during remediation. In particular, if a subset of particle types controls the progress of mass removal, both site characterization and performance modeling should be redirected accordingly. Future work will focus on clarifying the types of realistic conditions for which intraparticle diffusion and nonlinear sorption would be expected to control aquifer restoration in realistic multidimensional domains.

Extending the conceptual model and its numerical implementation to field scale applications will need to address several limitations. In particular, while the relevance of a detailed description of grain scale processes for field scale transport and remediation is an open question, the estimation and proper scale-up of isotherm and intraparticle diffusion parameters will be challenging. Furthermore, the computational burden of simulating intraparticle diffusion is likely to limit the scale of potential applications, even with the parallel-friendly algorithms presented in this work. However, there are several opportunities for further improvement to numerical performance, 32

including the use of an implicit solution technique for the intraparticle diffusion equations, adaptive time-stepping or other procedures to control operator-splitting errors, and/or more aggressive parallelization strategies (e.g., at the particle level).

33

APPENDIX A.1

Bisection algorithm

For problems characterized by nonlinear equilibrium sorption, an iterative solution is required to solve Eq. (10). Although Newton-Raphson iteration is an efficient approach for many problems of interest, convergence is not guaranteed, particularly for highly nonlinear sorption isotherms. As an alternative, a simple bisection procedure was implemented to identify the root to Eq. (10) located between pre-determined upper and lower bounds (Cupper and Clower, respectively). For the example problems presented in this paper, the initial bounds were established as follows:

(

C lower = min C * , C lT1

)

(

C upper = max C * , C lT1

(A.1a)

)

(A.1b)

where the C variables are as defined for Equation 10.

In a very few cases, the limits defined by Eq. A.1 will not contain the solution. In such cases, the following (wider) interval is recommended:

ρ   C lower = min  0, C * + b q T1  n  

(A.2a)

ρ   C upper = max  0, C * + b q T1  n  

(A.2b)

34

A.2

Finite difference formulation for diffusion equation

The finite difference formulations for Eqs. (11) are based on a forward time differencing and variable radial spacing scheme. The variable radial spacing scheme was designed to yield smaller radial spacing near the outer boundary. The non-uniform radial spacing was assigned such that the volumes of each spherical shell between the adjacent nodes were equal. A schematic representation of the node placement for spherical particles is shown in Figure A.1.

Figure A.1: Node placement scheme for a spherical particle.

For an unequally spaced discretization, the derivatives on the right hand side of Eq. (11) can be expressed for the kth particle at ith node with non-uniform radial spacing [55] as:

35

c − (1 + β k ,i )clk ,i + β k ,i clk ,i −1 ∂ 2 clk Dak = Dak lk ,i +1 2 ∂rk ∆rk ,i + ∆rk ,i − (1 + β k ,i ) 2

(A.3)

2 2 2 Dak ∂clk 2 Dak clk ,i +1 − (1 − β k ,i )clk ,i − β k ,i clk ,i −1 = rk ∂rk rk ,i (1 + β k ,i )∆rk ,i +

where, ∆rk ,i + = rk ,i +1 − rk ,i ; ∆rk ,i − = rk ,i − rk ,i −1 ; and β k ,i =

(A.4)

∆rk ,i+ . ∆rk ,i −

For the time step denoted by the superscript j, Eq. (11) reduces to:

clkj +,i1 − clkj ,i ∆t d

= Dak (bk ,i ,i +1clk ,i +1 + bk ,i ,i clk ,i + bk ,i ,i −1clk ,i −1 ) 2 Dak + (a k ,i ,i+1clk ,i+1 + a k ,i,i clk ,i + a k ,i ,i−1clk ,i−1 ) rk ,i

(A.5)

where, the coefficients of clk,i-1 , clk,i and clk,i+1 are:

bk ,i ,i +1 =

bk ,i ,i

1 ∆rk ,i + ∆rk ,i − (1 + β k ,i )

ak ,i ,i +1 =

1 ( 1 + β k ,i )∆rk ,i +

2

− (1 + β k ,i ) = ∆rk ,i+ ∆rk ,i− (1 + β k ,i )

ak ,i ,i

− ( 1 − β k2,i ) = ( 1 + β k ,i )∆rk ,i +

2 bk ,i ,i −1 =

β k ,i

a k ,i ,i −1 =

∆rk ,i + ∆rk ,i − (1 + β k ,i )

− β k2,i ( 1 + β k ,i )∆rk ,i +

(A.6)

2

Using equation A.5, the time derivatives at the inner nodes 2…(N-1) can be evaluated using the concentrations at the adjacent nodes.

36

At the center of particle, a Neumann boundary condition is applied at all the times; i.e. the diffusive flux at the center of particle is zero, yielding [56]:  ∂clk ,1  6 Dak  ∂t  = (∆r )2 (clk , 2 − clk ,1 ), at rk = 0, t > 0   k ,1+

(A.7)

In this work, the concentration at the particle boundary was developed by enforcing a mass balance on the elemental volume, with identical boundary concentration for all particles at exterior interface. The mass balance for a well-mixed system of elemental volume expressed as a function of boundary concentration for a reaction time step between time T2 and T2+1 can be expressed as:

(

nC

T2 +1 l

)

 rk ,i +1  T +1 m kT2,i++11 − m kT2,i+1   + ∑ ∑  N k ∫  m k 2,i + ( rk − rk ,i ) 4πrk2 drk  rk ,i +1 − rk ,i k i =1   rk ,i   

( )

= n C lT2

N −1

 rk ,i +1  T   m T2 − m kT2,i + ∑ ∑  N k ∫  m k 2,i + k ,i +1 ( rk − rk ,i ) 4π rk2 drk  rk ,i +1 − rk ,i k i =1   rk ,i   

(A.8)

N −1

T T T where, mk 2,i = n pk clk2,i + ρ bk qsk2,i , represents the “node mass” for particle k and node i at time T2; Nk

is the number of particles of identity k. The concentrations and nodal masses at time T2,

mkT2,i ∀ i = 1..N , are known from the previous time step and the concentrations/node masses at T +1 time T2+1, mk 2,i ∀ i = 1..(N − 1) , are known by solving Eq. A.5 and Eq. A.7.

Once the intraparticle diffusion equation is solved for the interior nodes, Equation A.8 can be solved for the only unknown which is concentration at exterior node Cl or clk,N at time T2+1 which is identical for all particles. Alternative formulations for the exterior boundary that use flux at boundary node for mass balance are discussed by Singh [57]. 37

NOMENCLATURE

Greek Letters ∆t

time interval between T1 and T1+1 [T]

∆td

time interval between T2 and T2+1 [T]

ρb

bulk density of the porous media [M/L3]

ρbk

bulk density for particle of identity k [M/L3]

∆r

intraparticle node spacing in diffusion domain [L]

∆x

grid spacing in transport domain [L]

Roman Letters A

Polanyi isotherm parameter

B

Polanyi isotherm parameter

b

Langmuir isotherm parameter

Cl

volume averaged aqueous phase solute concentration in [M/L3].

Cs

volume averaged sorbed solute concentration [M/L3].

C0

aqueous phase solute concentration at domain entrance [M/L3]

Cr

Courant number

cl

intraparticle aqueous phase solute concentration [M/L3]

D

dispersion coefficient [L2/T]

Da

effective intraparticle diffusion coefficient [L2/T]

Da

Damkohler number for intraparticle diffusion

F

mesh Fourier number

KP

partition coefficient in the Polanyi-partition isotherm [L3/M]

Kf

Freundlich isotherm parameter

nf

Freundlich isotherm parameter

n

aquifer porosity

np

intraparticle porosity

Pe

mesh Peclet number

Q0

Langmuir isotherm parameter

Qp

Polanyi isotherm parameter 38

qs

intraparticle sorbed mass fraction [M/M]

rk

radial distance from center of the particle of class k [L]

t

time [T]

v

groundwater velocity [L/T]

x

distance from the domain entrance [L]

Sw

solute aqueous solubility [M/L3]

Subscripts/Superscripts l

aqueous phase

s

sorbed phase

i

node number for discretized transport domain (x)

j

node number for discretized particles (r)

k

particle identity for particle with identical diffusion and sorption characteristics

T1

time sequence number for transport equation

T2

time sequence number for diffusion equation

Abbreviations ADRE

advective-dispersive-reactive-equation

HOC

hydrophobic organic contaminant

ODE

ordinary differential equation

MB

mass balance

PP

Polanyi-Partition

TCE

trichloroethylene

39

ACKNOWLEDGEMENTS

This study was funded by Strategic Environmental Research and Development Program (SERDP) project # ER-1738. Anshuman Singh’s doctoral studies were partially supported by a National Science Foundation IGERT traineeship award, number DGE-0333417. The computational support for this study was provided by the Center for Computational Research, University at Buffalo.

References [1] McCarty PL, M Reinhard, BE Rittmann. Trace organics in groundwater. Environmental Science & Technology. 15 (1981) 40-51, doi: 10.1021/es00083a003. [2] Roberts PV, PL McCarty, M Reinhard, J Schreiner. Organic contaminant behavior during groundwater recharge. Journal (Water Pollution Control Federation). 52 (1980) 161-72, doi: 10.2307/25040559. [3] Schwarzenbach RP, J Westall. Transport of nonpolar organic compounds from surface water to groundwater. Laboratory sorption studies. Environmental Science & Technology. 15 (1981) 1360-7, doi: 10.1021/es00093a009. [4] Legrand HE. Patterns of contaminated zones of water in the ground. Water Resour Res. 1 (1965) 83-95, doi: 10.1029/WR001i001p00083. [5] Chen W, AT Kan, CJ Newell, E Moore, MB Tomson. More Realistic Soil Cleanup Standards with Dual-Equilibrium Desorption. Ground Water. 40 (2002) 153-64, doi: 10.1111/j.17456584.2002.tb02500.x. [6] Sale T, C Newell. Guide for selecting remedies for subsurface releases of chlorinated solvents. Environmental Security Technology Certification Program (ESTCP Project ER200530) United States, 2011. pp. 146p. [7] Berglund S. The effect of Langmuir sorption on pump-and-treat remediation of a stratified aquifer. Journal of Contaminant Hydrology. 18 (1995) 199-220, doi: 10.1016/01697722(94)00057-O. [8] Berglund S, V Cvetkovic. Contaminant displacement in aquifers: Coupled effects of flow heterogeneity and nonlinear sorption. Water Resour Res. 32 (1996) 23-32, doi: 10.1029/95wr02767. 40

[9] Rabideau AJ, CT Miller. Two-dimensional modeling of aquifer remediation influenced by sorption nonequilibrium and hydraulic conductivity heterogeneity. Water Resour Res. 30 (1994) 1457-70, doi: 10.1029/93wr03414. [10] Rivett MO, SW Chapman, RM Allen-King, S Feenstra, JA Cherry. Pump-and-treat remediation of chlorinated solvent contamination at a controlled field-experiment site. Environmental Science & Technology. 40 (2006) 6770-81, doi: 10.1021/es0602748. [11] Weber Jr WJ, PM McGinley, LE Katz. Sorption phenomena in subsurface systems: Concepts, models and effects on contaminant fate and transport. Water Research. 25 (1991) 499528, doi: 10.1016/0043-1354(91)90125-A. [12] Mackay DM, JA Cherry. Groundwater contamination: pump-and-treat remediation. Environmental Science & Technology. 23 (1989) 630-6, doi: 10.1021/es00064a001. [13] Crittenden JC, NJ Hutzler, DG Geyer, JL Oravitz, G Friedman. Transport of organic compounds with saturated groundwater flow: Model development and parameter sensitivity. Water Resour Res. 22 (1986) 271-84, doi: 10.1029/WR022i003p00271. [14] Ball WP, PV Roberts. Long-term sorption of halogenated organic chemicals by aquifer material. 2. Intraparticle diffusion. Environmental Science & Technology. 25 (1991) 1237-49, doi: 10.1021/es00019a003. [15] Ball WP, PV Roberts. Long-term sorption of halogenated organic chemicals by aquifer material. 1. Equilibrium. Environmental Science & Technology. 25 (1991) 1223-37, doi: 10.1021/es00019a002. [16] Bahr JM. Analysis of nonequilibrium desorption of volatile organics during field test of aquifer decontamination. Journal of Contaminant Hydrology. 4 (1989) 205-22, doi: http://dx.doi.org/10.1016/0169-7722(89)90009-0. [17] Zhang Z, ML Brusseau. Nonideal transport of reactive solutes in heterogeneous porous media: 5. Simulating regional-scale behavior of a trichloroethene plume during pump-and-treat remediation. Water Resources Research. 35 (1999) 2921-35, doi: 10.1029/1999WR900162. [18] Zheng C, PP Wang. MT3DMS: A modular three-dimesional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; documentaion and user's guide. Contract Report SERDP-99-1, US Army Engineer Research and Development. US Army Corps of Engineers1999. [19] Huang W, TM Young, MA Schlautman, H Yu, WJ Weber. A distributed reactivity model for sorption by soils and sediments. 9. General isotherm nonlinearity and applicability of the dual reactive domain model. Environmental Science & Technology. 31 (1997) 1703-10, doi: 10.1021/es960677f.

41

[20] Xia G, WP Ball. Adsorption-partitioning uptake of nine low-polarity organic chemicals on a natural sorbent. Environmental Science & Technology. 33 (1998) 262-9, doi: 10.1021/es980581g. [21] Allen-King RM, P Grathwohl, WP Ball. New modeling paradigms for the sorption of hydrophobic organic chemicals to heterogeneous carbonaceous matter in soils, sediments, and rocks. Advances in Water Resources. 25 (2002) 985-1016, doi: 10.1016/S0309-1708(02)000453. [22] Freundlich HMF. Über die Adsorption in Löesungen. Zeitschrift für Physikalische Chemie. (1906) 385-470. [23] Langmuir I. The adsorption of gases on plane surfaces of glass, mica and platinum. Journal of the American Chemical Society. 40 (1918) 1361-403, doi: 10.1021/ja02242a004. [24] Clement TP. RT3D: A modular computer code for simulating reactive multispecies transport in 3-dimensional groundwater systems. Version 1.0 ed, Pacific Northwest National Laboratory, Richland, Washington, 1997. [25] Matott LS, AJ Rabideau. NIGHTHAWK-A program for modeling saturated batch and column experiments incorporating equilibrium and kinetic biogeochemistry. Computers & Geosciences. 36 (2010) 253-6, doi: 10.1016/j.cageo.2009.05.004. [26] Parkhurst DL, CAJ Appelo. User's guide to {PHREEQC} (version 2) - A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. U.S. Geological Survey1999. [27] Cheng H-P, G-T Yeh. Development and demonstrative application of a 3-D numerical model of subsurface flow, heat transfer, and reactive chemical transport: 3DHYDROGEOCHEM. Journal of Contaminant Hydrology. 34 (1998) 47-83, doi: 10.1016/S0169-7722(98)00084-9. [28] van Genuchten MT, FN Dalton. Models for simulating salt movement in aggregated field soils. Geoderma. 38 (1986) 165-83, doi: 10.1016/0016-7061(86)90013-3. [29] Pignatello JJ, B Xing. Mechanisms of slow sorption of organic chemicals to natural particles. Environmental Science & Technology. 30 (1995) 1-11, doi: 10.1021/es940683g. [30] Cunningham JA, CJ Werth, M Reinhard, PV Roberts. Effects of grain-scale mass transfer on the transport of volatile organics through sediments: 1. Model development. Water Resour Res. 33 (1997) 2713-26, doi: 10.1029/97WR02425. [31] Grathwohl P. Diffusion in natural porous media. 1 ed. Kluwer Academic Publishers, 1998.

42

[32] Miller CT, JA Pedit. Use of a reactive surface-diffusion model to describe apparent sorption-desorption hysteresis and abiotic degradation of lindane in a subsurface material. Environmental Science & Technology. 26 (1992) 1417-27, doi: 10.1021/es00031a021. [33] Kleineidam S, H Rügner, P Grathwohl. Desorption kinetics of phenanthrene in aquifer material lacks hysteresis. Environmental Science & Technology. 38 (2004) 4169-75, doi: 10.1021/es034846p. [34] Miller CT, WJ Weber Jr. Modeling the sorption of hydrophobic contaminants by aquifer materials—II. Column reactor systems. Water Research. 22 (1988) 465-74, doi: 10.1016/00431354(88)90041-3. [35] Werth CJ, JA Cunningham, PV Roberts, M Reinhard. Effects of grain-scale mass transfer on the transport of volatile organics through sediments: 2. Column results. Water Resour Res. 33 (1997) 2727-40, doi: 10.1029/97WR02426. [36] Xia G, WP Ball. Polanyi-Based Models for the Competitive Sorption of Low-Polarity Organic Contaminants on a Natural Sorbent. Environmental Science & Technology. 34 (2000) 1246-53, doi: 10.1021/es9812453. [37] Yeh GT, VS Tripathi. A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resour Res. 25 (1989) 93-108, doi: 10.1029/WR025i001p00093. [38] Rabideau AJ. MOUSER Version 1: Moderately User-friendly Reactive Transport Model (available from www.groundwater.buffalo.edu). Department of Civil, Structural, and Environmental Engineering, University at Buffalo, Buffalo, NY2003. [39] Pedit JA, CT Miller. Heterogeneous sorption processes in subsurface systems. 2. Diffusion modeling approaches. Environmental Science & Technology. 29 (1995) 1766-72, doi: 10.1021/es00007a012. [40] Liedl R, T Ptak. Modelling of diffusion-limited retardation of contaminants in hydraulically and lithologically nonuniform media. Journal of Contaminant Hydrology. 66 (2003) 239-59, doi: 10.1016/S0169-7722(03)00028-7. [41] Rasmuson A, I Neretnieks. Exact solution of a model for diffusion in particles and longitudinal dispersion in packed beds. AIChE Journal. 26 (1980) 686-90, doi: 10.1002/aic.690260425. [42] Rosen JB. Kinetics of a fixed bed system for solid diffusion into spherical particles. The Journal of Chemical Physics. 20 (1952) 387-94, doi: 10.1063/1.1700431. [43] Kanney JF, CT Miller, DA Barry. Comparison of fully coupled approaches for approximating nonlinear transport and reaction problems. Advances in Water Resources. 26 (2003) 353-72, doi: 10.1016/S0309-1708(02)00188-4. 43

[44] Miller CT, AJ Rabideau. Development of split-operator, Petrov-Galerkin methods to simulate transport and diffusion problems. Water Resour Res. 29 (1993) 2227-40, doi: 10.1029/93wr00528. [45] Barry DA, K Bajracharya, CT Miller. Alternative split-operator approach for solving chemical reaction/groundwater transport models. Advances in Water Resources. 19 (1996) 26175, doi: 10.1016/0309-1708(96)00002-4. [46] Barry DA, K Bajracharya, M Crapper, H Prommer, CJ Cunningham. Comparison of splitoperator methods for solving coupled chemical non-equilibrium reaction/groundwater transport models. Mathematics and Computers in Simulation. 53 (2000) 113-27, doi: 10.1016/S03784754(00)00182-8. [47] Valocchi AJ, M Malmstead. Accuracy of operator splitting for advection-dispersionreaction problems. Water Resour Res. 28 (1992) 1471-6, doi: 10.1029/92wr00423. [48] Steefel CI, SB Yabusaki. OS3D/GIMRT: Software for Modeling MulticomponentMultidimensional Reactive Transport. Pacific Noethwest National Laboratory, Washington, USA1996. [49] Gasda, S. E., MW Farthing, CE Kees, CT Miller. Adaptive split operator methods for modeling transport phenomena in porous media syetems. Advances in Water Resources. 34 (2011) 1268-1282. [50] Zheng C, GD Bennett. Applied contaminant transport modeling. 2nd Edition ed. WileyInterscience, 2002. [51] Grove DB, KG Stollenwerk. Computer model of one-dimensional equilibrium controlled sorption process. U.S. Geological Survey1984. [52] Munger ZW. Modeling nonlinear sorption of trichloroethene in natural sorbents with kerogen [MS]: University at Buffalo, State University of New York; 2012. [53] Matott LS, AJ Rabideau. ISOFIT - A program for fitting sorption isotherms to experimental data. Environmental Modelling & Software. 23 (2008) 670-6, doi: 10.1016/j.envsoft.2007.08.005. [54] Barry DA, CT Miller, PJ Culligan, K Bajracharya. Analysis of split operator methods for nonlinear and multispecies groundwater chemical transport models. Mathematics and Computers in Simulation. 43 (1997) 331-41. [55] Hoffman JD. Numerical methods for engineers and scientists. 2nd ed. CRC Press, 2001. [56] Crank J. Mathematics of Diffusion. 2 ed. Oxford Science Publications, 1975.

44

[57] Singh A. The effect of nonideal sorption in vapor and groundwater transport: Statistical and modeling tools. Doctoral dissertation: University at Buffalo, State University of New York; 2012.

45

Table 1: Simulation parameters for Figures 2 – 6 Transport Parameters v D n ρb Source

0.091 m/day 7.3284e-03m2/day 0.33 1.792 g/ml3 30 µg/L pulse for 10 days

TCE Isotherm Parameters [51] Kf = 4.371 µg/kg (L/µg)nf nf = 0.793 Q0= 10815 µg/kg b= 1.443E-04 µg/L Qp= 36561.4 µg/kg A= -0.318 B = 1.418 Sw = 1.4e+06 µg/L KP = 0.221 L/µg

Freundlich Langmuir Polanyi-partitioning

Discretization No. of spatial nodes

1001

No. of time steps for Courant number 0.05, 0.1, 0.5, 1.0 and total simulation time of 600 days

24000, 12000, 2400, 1200

Table 2: Additional parameters for intraparticle diffusion simulations (Figures 7 – 9) Radius,Rk Intraparticle porosity, npk Effective diffusion coefficient, Dak Bulk density of particles, ρbk

(i) 0.7 mm; (ii) 1.5 mm; (iii) 2.5 mm 0.0603 1.11e-06 m2/day 2.56 g/cm3

46

Table 3: Spatial and temporal discretization for intraparticle diffusion simulations (Figures 7 – 9) Isotherm Freundlich Polanyipartitioning

Number of particle nodes

Number of diffusion time steps

Particle Radius

Cr=0.1

Cr=0.5

Cr=1

Cr=0.1

Cr=0.5

Cr=1

0.7 mm 1.5 mm 2.5 mm 0.7 mm 1.5 mm 2.5 mm

28 45 14 28 46 14

13 21 7 14 22 7

10 15 6 10 16 6

317 306 328 317 342 328

303 310 327 357 342 327

334 299 442 334 344 442

Table 4: Simulation parameters multi-particle scenarios (Figures 10 – 11) Transport Parameters v D n Domain length Initial aqueous concentration

1.0 m/day 4.0e-03 m2/day 0.4 0.2 m 40 µg/L

Particle Parameters Sandstone

Dark Limestone Mass fraction 0.656 0.016 Radius (cm) 0.015 0.015 Number of intraparticle nodes 16 101 and 21* Intraparticle porosity 0.0195 0.0035 3 Bulk density (g/cm ) 2.63 2.73 2 Diffusion coefficient (m /d) 4.12e-06 8.38e-8 and 2.47e-06* nf Freundlich Kf (µg/kg (L/µg) ) 30.36 53.75 3683.85 Freundlich nf 0.66 0.67 0.33 * Two variations of the diffusion coefficients and discretization nodes were used for dark limestone.

47

Light Limestone 0.328 0.015 21 0.0054 2.715 2.47e-06

Highlights

•

We present split-operator models for solute transport with nonideal sorption.

•

Numerical experiments address transport with “dual mode” and nonlinear isotherms.

•

Numerical experiments address intraparticle diffusion for multi-particle systems.

48