An analytical model for simulating two-dimensional multispecies plume migration

The two-dimensional advection-dispersion equations coupled with sequential first-order decay reactions involving arbitrary number of species in groundwater system is considered to predict the two-dimensional plume behavior of decaying contaminant such as radionuclide and dissolved chlorinated solvent. Generalized analytical solutions in compact format are derived through the sequential application of the Laplace, finite Fourier cosine, and generalized integral transform to reduce the coupled partial differential equation system to a set of linear algebraic equations. The system of algebraic equations is next solved for each species in the transformed domain, and the solutions in the original domain are then obtained through consecutive integral transform inversions. Explicit form solutions for a special case are derived using the generalized analytical solutions and are compared with the numerical solutions. The analytical results indicate that the analytical solutions are robust, accurate and useful for simulation or screening tools to assess plume behaviors of decaying contaminants.


Introduction
Experimental and theoretical studies have been undertaken to understand the fate and transport of dissolved hazardous substances in subsurface environments because human health is threatened by a wide spectrum of contaminants in groundwater and soil.Analytical models are essential and efficient tools for understanding pollutants behavior in subsurface environments.Several analytical solutions for single-species transport problems have been reported for simulating the transport of various contaminants (Batu, 1989(Batu, , 1993(Batu, , 1996;;Chen et al., 2008aChen et al., , b, 2011;;Gao et al., 2010Gao et al., , 2012Gao et al., , 2013;;Leij et al., 1991Leij et al., , 1993;;Park and Zhan, 2001;Pérez Guerrero and Skaggs, 2010;Pérez Guerrero et al., 2013;van Genuchten and Alves, 1982;Yeh, 1981;Zhan et al., 2009;Ziskind et al., 2011).Transport processes of some contaminants such as radionuclides, dissolved chlorinated solvents and nitrogen generally involve a series of first-order or pseudo first-order sequential decay chain reactions.During migrations of decaying contaminants, mobile and toxic successor products may sequentially form and move downstream with elevated concentrations.Single-species analytical models do not permit transport behaviors of successor species of these decaying contaminants to be evaluated.Analytical models for multispecies transport equations coupled with first-order sequential decay reactions are useful tools for synchronous determination of the fate and transport of the predecessor and successor species of decaying contaminants.However, there are few analytical solutions for coupled multispecies transport equations compared to a large body of analytical solutions in the literature pertaining to the single-species advectivedispersive transport subject to a wide spectrum of initial and boundary conditions.
Mathematical approaches have been proposed in the literature to derive a limited number of one-dimensional analytical solutions or semi-analytical solutions for multispecies advective-dispersive transport equations sequentially coupled with first-order decay reactions.These include direct integral transforms with sequential substitutions (Cho, Published by Copernicus Publications on behalf of the European Geosciences Union. 1971; Lunn et al., 1996;van Genuchten, 1985;Mieles and Zhan, 2012), decomposition by change-of-variables with the help of existing single-species analytical solutions (Sun and Clement, 1999;Sun et al., 1999a, b), Laplace transform combined with decomposition of matrix diagonization (Quezada et al., 2004;Srinivasan and Clememt, 2008a, b), decomposition by change-of-variables coupled with generalized integral transform (Pérez Guerrero et al., 2009, 2010), sequential integral transforms in association with algebraic decomposition (Chen et al., 2012a, b).
Multi-dimensional solutions are needed for real world applications, making them more attractive than onedimensional solutions.Bauer et al. (2001) presented the first set of semi-analytical solutions for one-, two-, and three-dimensional coupled multispecies transport problem with distinct retardation coefficients.Explicit analytical solutions were derived by Montas (2003) for multi-dimensional advective-dispersive transport coupled with first-order reactions for a three-species transport system with distinct retardation coefficients of species.Quezada et al. (2004) extended the Clement (2001) strategy to obtain Laplace-domain solutions for an arbitrary decay chain length.Most recently, Sudicky et al. (2013) presented a set of semi-analytical solutions to simulate the three-dimensional multi-species transport subject to first-order chain-decay reactions involving up to seven species and four decay levels.Basically, their solutions were obtained species by species using recursion relations between target species and its predecessor species.For a straight decay chain, they derived solutions for up to four species and no generalized expressions with compact formats for any target species were obtained.Note that their solutions were derived for the first-type (Dirichlet) inlet conditions which generally bring about physically improper mass conservation and significant errors in predicting the concentration distributions especially for a transport system with a large longitudinal dispersion coefficient (Barry and Sposito, 1988;Parlange et al., 1992).Moreover, in addition to some special cases, the numerical Laplace transforms are required to obtain the original time domain solution.Besides the straight decay chain, the analytical model by Clement (2001) and Sudicky et al. (2013) can account for more complicated decay chain problems such as diverging, converging and branched decay chains.
Based on the aforementioned reviews, this study presents a parsimonious explicit analytical model for two-dimensional multispecies transport coupled by a series of first-order decay reactions involving an arbitrary number of species in groundwater system.The derived analytical solutions have four salient features.First, the third-type (Robin) inlet boundary conditions which satisfy mass conservation are considered.Second, the solution is explicit, thus solution can be easily evaluated without invoking the numerical Laplace inversion.Third, the generalized solutions with parsimonious mathematical structures are obtained and valid for any species of a decay chain.The parsimonious mathematical structures of the generalized solutions are easy to code into a computer program for implementing the solution computations for arbitrary target species.Fourth, the derived solutions can account for any decay chain length.The explicit analytical solutions have applications for evaluation of concentration distribution of arbitrary target species of the real-world decaying contaminants.The developed parsimonious model is robustly verified with three example problems and applied to simulate the multispecies plume migration of dissolved radionuclides and chlorinated solvent.

Derivation of analytical solutions
This study considers the problem of decaying contaminant plume migration.The source zone is located in the upstream of groundwater flow.The source zone can represent leaching of radionuclide from a radioactive waste disposal facility or release of chlorinated solvent from the residual NAPL phase into the aqueous phase.After these decaying contaminants enter the aqueous phase, they migrate by one-dimensional advection with flowing groundwater and by simultaneously longitudinal and transverse dispersion processes.While migrating in the groundwater system, the contaminants undergo linear isothermal equilibrium sorption and a series of sequential first-order decaying reactions.Sudicky et al. (2013) provided the detailed modeling scenario.The scenario considered in this study can be ideally described as shown in Fig. 1.A steady and uniform velocity in the x direction is considered in Fig. 1.The governing equations describing two-dimensional reactive transport of the decaying contaminants and their successor species undergoing linear isothermal equilibrium sorption and a series of sequential first-order decaying reactions can be mathematically written as where C i (x, y, t) is the aqueous concentration of species i is the retardation coefficient of species i [−].Note that these equations consider that the decay reactions occur simultaneously in both the aqueous and sorbed phases.If the decay reactions occur only in the aqueous phase, the retardation coefficients in the decay terms in the right-hand sides of Eqs. ( 1a) and (1b) become unity.For such a case, k i and k i−1 in the left-hand sides could be modified as k i R i and k i−1 R i−1 to facilitate the application of the derived analytical solutions obtained by Eqs.(1a) and (1b).
The initial and boundary conditions for solving Eqs.(1a) and (1b) are where f i (t) is the arbitrary time-dependent source concentration of species i applied at the source segment (H (y − y 1 )−H (y −y 2 )) at boundary (x = 0) which will be specified later [L], H (•) is the Heaviside function, L and W are the length and width of the transport system under consideration [L].Equation (2) implies that the transport system is free of solute mass at the initial time.Equation (3) means that a third-type boundary condition satisfying mass conservation at the inlet boundary is considered.Equation (4) considers the concentration gradient to be zero at the exit boundary based on the mass conservation principle.Such a boundary condition has been widely  used for simulating solute transport in a finite-length system.Equations ( 5) and (6) assume no solute flux across the lower and upper boundaries.It is noted that in Eq. (3), we assume arbitrary time-dependent sources of species i uniformly distributed at the segment (y 1 ≤ y ≤ y 2 ) of the inlet boundary (x = 0), the so-called Heaviside function source concentration profile.Relative to the first type boundary conditions used by Sudicky et al. (2013), the third-type boundary conditions which satisfy mass conservation at the inlet bound- ary (Barry and Sposito, 1988;Parlange et al., 1992) are used herein.Sudicky et al. (2013) considered the source concentration profiles as Gaussian or Heaviside step functions.If Gaussian distributions are desired, we can easily replace the Heaviside function in the right-hand side of Eq. (3) with a Gaussian distribution.
Equations ( 1)-( 6) can be expressed in dimensionless form as where Our solution strategy used is extended from the approach proposed by Chen at al. (2012a, b).The core of this approach is that the coupled partial differential equations are converted into an algebraic equation system via a series of integral transforms and the solutions in the transformed domain for each species are directly and algebraically obtained by sequential substitutions.
Following Chen et al. (2012a, b), the generalized analytical solutions in compact formats can be obtained as follows (with detailed derivation provided in Appendix A)   (1985).25044  230 Th, i = 3 0.443684 × 10 −3  0.593431  −0.593874  226 Ra, i = 4 −0.516740 × 10 −6 0.120853 × 10 − and where Concise expressions for arbitrary target species such as described in Eqs. ( 13) to (15) facilitate the development of a computer code for implementing the computations of the analytical solutions.
The generalized solutions of Eq. ( 13) accompanied by two corresponding auxiliary functions p i (ξ l , n, T ) and q i (ξ l , n, T )in Eqs. ( 14)-( 15) can be applied to derive analytical solutions for some special-case inlet boundary sources.
Here the time-dependent decaying source which represents the specific release mechanism defined by the Bateman equations (van Genuchten, 1985) is considered.A Bateman-type source is described by or in dimensionless form, The coefficients b im and δ m = µ m + γ m account for the firstorder decay reaction rate (µ m ) of each species in the waste source and the release rate (γ m ) of each species from the waste source, λ m = δ m L v .By substituting Eq. (16b) into Eqs.( 13)-( 15), we obtain where and

Convergence behavior of the Bateman-type source solution
Based on the special-case analytical solutions in Eq. ( 17) supported by two auxiliary functions, defined in Eqs. ( 18) and ( 19), a computer code was developed in FORTRAN 90 language with double precision.The details of the FORTRAN computer code are described in Supplement.
The derived analytical solutions in Eqs. ( 17)-( 19) consist of summations of double infinite series expansions for the finite Fourier cosine and generalized integral transform inversions, respectively.It is straightforward to sum up these two infinite series expansions term by term.To avoid time-consuming summations of these infinite series expansions, the convergence tests should be routinely executed to determine the optimal number of the required   terms for evaluating analytical solutions to the desired accuracies.Two-dimensional four-member radionuclide decay chain 238 Pu → 234 U → 230 Th → 226 Ra is considered herein as convergence test example 1 to demonstrate the convergence behavior of the series expansions.This convergence test example 1 is modified from a one-dimensional radionuclide decay chain problem originated by Higashi and Pigford (1980) and later applied by van Genuchten (1985) to illustrate the applicability of their derived solution.The important model parameters related to this test example are listed in Tables 1 and 2. The inlet source is chosen to be symmetrical with respect to the x-axis and conveniently arranged in the 40 m ≤ y ≤ 60 m segment at the inlet boundary.
In order to determine the optimal term number of series expansions for the finite Fourier cosine transform inverse to achieve accurate numerical evaluation, we specify a sufficiently large number of series expansions for the generalized transform inverse so that the influence of the number of series expansions for the generalized integral transform inverse on convergence of series expansion for finite Fourier cosine transform inverse can be excluded.A similar concept is used when investigating the required number of terms in the series expansions for the generalized integral transform inverse.An alternative approach is conducted by simultaneously varying the term numbers of series expansions for the generalized integral transform inverse and the finite Fourier cosine transform inverse.
Tables 3, 4 and 5 give results of the convergence tests up to 3 decimal digits of the solution computations along the three transects (inlet boundary at x = 0 m, x = 25 m, and exit boundary at x = 250 m).In these tables M and N are defined as the numbers of terms summed for the generalized integral transform inverse and finite Fourier cosine transform inverse, respectively.It is observed that M and N are related closely to the true values of the solutions.For smaller true values, the solutions must be computed with greater M and N .However, convergences can be drastically speeded up if lower calculation precision (e.g. 2 decimal digits accuracy) is acceptable.For example, (M, N ) = (100 200) is sufficient for 2 decimal digits accuracy, while for 3 decimal digits accuracy we need (M, N ) = (1600, 8000).Two decimal digits accuracy is acceptable for most practical problems.It is also found that M increases and N decreases with increasing x.
To further examine the series convergence behavior, example 2 considers a transport system of large aspect ratio ( L W = 2500 m 100 m ) and a narrower source segment, 45 m ≤ y ≤ 55 m, on the inlet boundary.Tables 6 and 7 present results of the convergence tests of the solution computations along two transects (inlet boundary and x = 250 m).Tables 6 and 7 also show similar results for the dependences of M and N on x.
Note that larger M and N are required for each species in this test example, suggesting that the evaluation of the solution for a large aspect ratio requires more series expansion terms to achieve the same accuracy as compared to example 1. Detailed results of the convergence test examples 1 and 2 are provided in the Supplement.
Using the required numbers determined from the convergence test, the computational time for evaluation of the solutions at 50 different observations only takes 3.782, 11.325, 23.95 and 67.23 s computer clock time on an Intel Core i7-2600 3.40 MHz PC for species 1, 2, 3, and 4 in the comparison of example 1.

Comparison of the analytical solutions with the numerical solutions
Three comparison examples are considered to examine the correctness and robustness of the analytical solutions and the      Moridis and Reddell (1991).A computer code for the LTFD solution is written in FORTTRAN language with double precision.The details of the FORTRAN computer code are described in Supplement.Figures 2, 3 and 4 depicts the spatial concentration distribution along one longitudinal direction (y = 50 m) and two transverse directions (x = 0 m and x = 25 m) for convergence test example 1 at t = 1000 years obtained from analytical solutions and numerical solutions.Figures 5, 6 and 7 present the spatial concentration distribution along one longitudinal direction (y = 50 m) and two transverse directions (x = 0 m and x = 25 m) for the convergence test example 2 at t = 1000 years obtained from analytical solutions and numerical solutions.Excellent agreements between the two solutions for both examples are observed for a wide spectrum of concentration, thus warranting the accuracy and robustness of the developed analytical model.
The third example involves a 10 species decay chain previously presented by Srinivasan and Clement (2008a) to evaluate the performance of their one-dimensional analytical solutions.The relevant model parameters are summarized in Tables 8 and 9. Our computer code is also compared against the LTFD solutions for this example.Figure 8 depicts the spatial concentration distribution at t = 20 days obtained analytically and numerically.Again there is excellent agreement between the analytical and numerical solutions, demonstrating the performance of our computer code for simulating transport problems with a long decay chain.The three comparison results clearly establish the correctness of the analytical model and the accuracy and capability of the computer code.

Assessing physical and chemical parameters on the radionuclide plume migration
Physical processes and chemical reactions affect the extent of contaminant plumes, as well as concentration levels.To illustrate how the physical processes and chemical reac- Figure 9 illustrates the spatial concentration of four species at t = 1000 years for the three sets of dispersion coefficients.The mobility of plumes of 234 U and 230 Th is retarded because of their stronger sorption ability.Hence the least retarded 226 Ra plume extensively migrated to 200 m × 60 m area in the simulation domain, whereas the 234 U and 230 Th Table 6.Solution convergence of each species concentration at transect of inlet boundary (x = 0 m) for four-species radionuclide transport problem considering simulated domain of L = 2500 m, W = 100 m subject to Bateman-type sources located at 45 m ≤ y ≤ 55 m for t = 1000 years (M = number of terms summed for inverse generalized integral transform; N = number of terms summed for inverse finite Fourier cosine transform).When we investigate the required M for inverse generalized integral transform, N = 12 800 for the finite Fourier cosine transform inverse are used.When we investigate the required N for inverse finite Fourier cosine transform, M = 6400 for the generalized transform inverse are used. 238Pu  7. Solution convergence of each species concentration at transect of x = 250 m for four-species radionuclide transport problem considering simulated domain of L = 2500 m, W = 100 m subject to Bateman-type sources located at 45 m ≤ y ≤ 55 m for t = 1000 years (M = number of terms summed for inverse generalized integral transform; N = number of terms summed for inverse finite Fourier cosine transform).When we investigate the required M for inverse generalized integral transform, N = 160 for the finite Fourier cosine transform inverse are used.When we investigate the required N for inverse finite Fourier cosine transform, M = 12 800 for the generalized transform inverse are used. 238Pu   , 2, . . . , 10 3, 2, 1.5, 1.25, 2.75, 1, 0.75, 0.5, 0.25, 0.1 Source decay constant, λ m [year −1 ] m = 1, 2,. . .,10 0.1, 0.75, 0.5, 0.25, 0, 0, 0.3, 1, 0, 0.65 Table 9. Coefficients of Bateman-type boundary source for ten-species transport problem used by Srinivasan and Clement (2008b).The moderate mobility of 238 Pu reflects the fact that it is a medial sorbed member of this radionuclide decay chain.The high concentration level of 234 U accounts for the high first-order decay rate constant of its parent species 238 Pu and its own low first-order decay rate constant.The plume extents and concentration levels may be sensitive to longitudinal and transverse dispersion.Increase of the longitudinal and/or transverse dispersion coefficients enhances the spreading of the plume extensively along the longitudinal and/or transverse directions, thereby lowering the plume concentration level.Because the concentration levels of the four radionuclides are influenced by both source release rates and decay chain reactions, 230 Th has the least extended plume area, while 226 Ra has the greatest plume area for all three sets of dispersion coefficients.These dispersion coefficients only affect the size of plumes of the four radionuclide, but the order of their relative plume size remains the same (i.e. 226Ra > 238 Pu > 234 U > 230 Th for the simulated condition).Indeed, in the reactive contaminant transport, the chemical parameters of sorption and decay rate are more important than the physical parameters of dispersion coefficients that govern the order of the plume extents and the concentration levels.

Simulating the natural attenuation of chlorinated solvent plume migration
Natural attenuation is the reduction in concentration and mass of the contaminant due to naturally occurring processes in the subsurface environment.The process is monitored for regulatory purposes to demonstrate continuing attenuation of the contaminant reaching the site-specific regulatory goals within reasonable time, hence, the use of the term monitored natural attenuation (MNA).MNA has been widely accepted as a suitable management option for chlorinated solvent contaminated groundwater.Mathematical models are widely used to evaluate the natural attenuation of plumes at chlorinated solvent sites.The multispecies transport analytical model developed in this study provides an effective tool for evaluating performance of the monitoring natural attenuation of plumes at a chlorinated solvent site because a series of daughter products produced during biodegradation of chlorinated solvent such as PCE → TCE → DCE → VC → ETH.Thus simulation of the natural attenuation of plumes a chlorinated solvent constitutes an attractive field application example of our multispecies transport model.
A study of 45 chlorinated solvent sites by McGuire et al. (2004) found that mathematical models were used at 60 % of these sites and that the public domain model BIOCHLOR  (Aziz et al., 2000) provided by the Center for Subsurface Modeling Support (CSMoS) of USEPA was the most commonly used model.An illustrated example from BIOCHLOR manual (Aziz et al., 2000) is considered to demonstrate the application of the developed analytical model.This example application demonstrated that BIOCHLOR can reproduce plume movement from 1965 to 1998 at the contaminated site of Cape Canaveral Air Station, Florida.The simulation conditions and transport parameters for this example application are summarized in Table 10.Constant source concentrations rather than exponentially declining source concentration of five-species chlorinated solvents are specified in the 90.7 m ≤ y ≤ 122.7 m segment at the inlet boundary (x = 0).This means that the exponents (λ im ) of Batemantype sources in Eqs.(16a) or (16b) need to be set to zero for the constant source concentrations and source intensity constants (b im ) are set to zero when subscript i does not equal to subscript m.Table 11 lists the coefficients of Bateman-type boundary source used for this example application involving the five-species dissolved chlorinated solvent problem.Spatial concentration contours of five-species at t = 1 year obtained from the derived analytical solutions for natural attenuation of chlorinated solvent plumes are depicted in Fig. 10.
It is observed that the mobility of plumes is quite sensitive to the species retardation factors, whereas the decay rate constants determine the plume concentration level.The plumes can migrate over a larger region for species having a low retardation factor such as VC.The low decay rate constants such as ETH have higher concentration distribution than the VC.It should be noted that a larger extent of plume observed for ETH in Fig. 10 is mainly attributed to the plume mass accumulation from the predecessor species VC that have a larger plume extent.The effect of high retardation of the ETH is hindered by the mass accumulation of the predecessor species VC.

Conclusions
We present an analytical model with a parsimonious mathematical format for two-dimensional multispecies advectivedispersive transport of decaying contaminants such as radionuclides, chlorinated solvents and nitrogen.The developed model is capable of accounting for the temporal and spatial development of an arbitrary number of sequential first-order decay reactions.The solution procedures involve applying a series of Laplace, finite Fourier cosine and generalized integral transforms to reduce a partial differential equation system to an algebraic system, solving for the algebraic system for each species, and then inversely transforming the concentration of each species in transformed domain into the original domain.Explicit special solutions for Bateman type source problems are derived via the generalized analytical solutions.The convergence of the series expansion of the generalized analytical solution is robust and accurate.These explicit solutions and the computer code are compared with the results computed by the numerical solutions.The two solutions agree well for a wide spectrum of concentration variations for three test examples.The analytical model is applied to assess the plume development of radionuclide and dissolved chlorinated solvent decay chain.
The results show that dispersion only moderately modifies the size of the plumes, without altering the relative order of the plume sizes of different contaminant.It is suggested that retardation coefficients, decay rate constants and the prede-cessor species plume distribution mainly govern the order of plume size in groundwater.Although there are a number of numerical reactive transport models that can account for multispecies advective-dispersive transport, our analytical model with a computer code that can directly evaluate the two-dimensional temporal-spatial concentration distribution of arbitrary target species without involving the computation of other species.The analytical model developed in this study effectively and accurately predicts the two-dimensional radionuclide and dissolved chlorinated plume migration.It is a useful tool for assessing the ecological and environmental impact of the accidental radionuclide releases such as the Fukushima nuclear disaster where multiple radionuclides leaked through the reactor, subsequently contaminating the local groundwater and ocean seawater in the vicinity of the nuclear plant.It is also a screening model that simulates remediation by natural attenuation of dissolved solvents at chlorinated solvent release sites.It should be noted the derived analytical model still has its application limitations.The model cannot handle the site with non-uniform groundwater flow or with multiple distinct zones.Furthermore, the developed model cannot simulate the more complicated decay chain problems such as diverging, converging and branched decay chains.The analytical model for more complicated decay chain problems can be pursued in the near future.
= −f i (T ) + β i e −α i T T 0 f i (τ )e α i τ dτ.(A31) The Laplace inverse of Q i (ξ l , n, s) can be also approached using the similar method.By taking Laplace inverse transform on Q i (ξ l , n, s), we have q i (ξ l , n, T ) = L −1 [Q i (ξ l , n, s)] as the summation of partial fractions and applying the inverse Laplace transform formula, one gets e −α i−j 1 T j 3 =i j 3 =i−k−1,j 3 =i−j 1 (α j 3 − α i−j 1 ) Recall that the inverse Laplace transform of F i−k−1 (s) is f i−k−1 (T ).Thus, the Laplace inverse transform of 1 j 2 =k+1 j 2 =0 s+α i−j 2 F i−k−1 (s) in Eq. ( 1) can be achieved using the convolution integral equation as e −α i−j 1 T T 0 e α i−j 1 τ f i−k−1 (τ )dτ j 3 =i j 3 =i−k−1,j 3 =i−j 2 α j 3 − α i−j 2 (A34) Putting Eq. (A34) into Eq.(A2) we can obtain the following form: e −α i−j 1 T T 0 e α i−j 1 τ f i−k−1 (τ )dτ j 3 =i j 3 =i−k−1,j 3 =i−j 2 α j 3 − α i−j 2 . (A35) Thus, the final solution can be expressed as Eq. ( 13) with the corresponding functions defined in Eqs. ( 14) and ( 15).Note that Eq. ( A33) is invalid for some of α i−j 2 being identical.For such conditions, we can still reduce 1 j 2 =k+1 j 2 =0 s+α i−j 2 to a sum of partial fraction expansion.However, it will lead to different Laplace inverse formulae.For example, the following formulae is used for all α i−j 2 being identical The generalized formulae for the cases with some of α i−j 2 being identical will not be provided herein because there are a large number of combinations of α i−j 2 .We suggest that the readers can pursue the solutions by following the similar steps for such specific conditions case by case.
The Supplement related to this article is available online at doi:10.5194/hess-20-733-2016-supplement.

Figure 1 .
Figure 1.Schematic representation of two-dimensional transport of decaying contaminants in a uniform flow field with flux boundary source located at of the inlet boundary.
Fig. 2.Figure 2. Comparison of spatial concentration profiles of four species along the longitudinal direction (= 50 m) at t = 1000 years obtained from derived analytical solutions and numerical solutions for convergence test example 1 of four-member radionuclide decay chain 238 Pu → 234 U → 230 Th → 226 Ra.

Figure 2 .
Fig. 2.Figure 2. Comparison of spatial concentration profiles of four species along the longitudinal direction (= 50 m) at t = 1000 years obtained from derived analytical solutions and numerical solutions for convergence test example 1 of four-member radionuclide decay chain 238 Pu → 234 U → 230 Th → 226 Ra.

Figure 3 .
Fig. 3. Figure 3.Comparison of spatial concentration profiles of four species along the transverse direction (= 0 m) at t = 1000 years obtained from derived analytical solutions and numerical solutions for convergence test example 1 of four-member radionuclide decay chain 238 Pu → 234 U → 230 Th → 226 Ra.

Figure 4 .
Figure 4. Comparison of spatial concentration profiles of four species along the transverse direction (= 25 m) at t = 1000 years obtained from derived analytical solutions and numerical solutions for convergence test example 1 of four-member radionuclide decay chain 238 Pu → 234 U → 230 Th → 226 Ra.

Figure 5 .
Fig. 5 Figure 5.Comparison of spatial concentration profiles of four species along the longitudinal direction (= 50 m) at t = 1000 years obtained from derived analytical solutions and numerical solutions for convergence test example 2 of four-member radionuclide decay chain 238 Pu → 234 U → 230 Th → 226 Ra.

Figure 6 .
Fig. 6 Figure 6.Comparison of spatial concentration profiles of four species along the transverse direction (= 0 m) at t = 1000 years obtained from derived analytical solutions and numerical solutions for convergence test example 2 of four-member radionuclide decay chain 238 Pu → 234 U → 230 Th → 226 Ra.
Fig. 7.Figure 7. Comparison of spatial concentration profiles of four species along the transverse direction (= 25 m) at t = 1000 years obtained from derived analytical solutions and numerical solutions for convergence test example 2 of four-member radionuclide decay chain 238 Pu → 234 U → 230 Th→ 226 Ra. e

Figure 8 .
Fig. 8Figure8.Comparison of spatial concentration profiles of ten-species along x direction at t = 20 days obtained from derived analytical solutions and numerical solutions for the test example 3 of ten species decay chain used bySrinivasan and Clement (2008b).

Figure 9 .
Figure 9. Effects of physical processes and chemical reactions on the concentration contours of four-species at t = 1000 years obtained from derived analytical solutions for four-member decay chain 238 Pu → 234 U → 230 Th→ 226 Ra.

Figure 10 .
Figure 10.Spatial concentration contours of five-species at t = 1 year obtained from derived analytical solutions for natural attenuation of chlorinated solvent plumes PCE → TCE → DCE → VC → ETH.
tions affect multispecies plume development, we consider the four-member radionuclide decay chain used in the previous convergence test and solution verification.The model parameters are the same, except that the longitudinal (D L ) and transverse (D T ) dispersion coefficients are varied.Three sets of longitudinal and transverse dispersion coefficients D L = 1000, D T = 100; D L = 1000, D T = 200; D L = 2000, D T = 200 (all in m 2 year −1 ) are tested, all for a simulation time of 1000 years.

Table 1 .
Transport parameters used for convergence test example 1 involving the four-species radionuclide decay chain problem used byvan Genuchten (1985).

Table 2 .
Values for coefficients of Bateman-type boundary source for four-species transport problem used by van Genuchten

Table 3 .
Solution convergence of each species concentration at transect of inlet boundary (x = 0) for four-species radionuclide transport problem considering simulated domain of L = 250 m, W = 100 m, subject to Bateman-type sources located at 40 m ≤ y ≤ 60 m for t = 1000 s year (M = number of terms summed for inverse generalized integral transform; N = number of terms summed for inverse finite Fourier cosine transform).When we investigate the required M for inverse generalized integral transform, N = 16 000 for the finite Fourier cosine transform inverse are used.When we investigate the required N for inverse finite Fourier cosine transform, M = 1600 for the generalized transform inverse are used.

Table 4 .
Solution convergence of each species concentration at transect of x = 25 m for four-species radionuclide transport problem considering simulated domain of L = 250 m, W = 100 m, subject to Bateman-type sources located at 40 m ≤ y ≤ 60 m for t = 1000 year (M = number of terms summed for inverse generalized integral transform; N = number of terms summed for inverse finite Fourier cosine transform).When we investigate the required M for inverse generalized integral transform, N = 160 for the finite Fourier cosine transform inverse are used.When we investigate the required N for inverse finite Fourier cosine transform, M = 1600 for the generalized transform inverse are used.

Table 5 .
Solution convergence of each species concentration at transect of exit boundary (x = 250 m) for four-species radionuclide transport problem considering simulated domain of L = 250 m, W = 100 m subject to Bateman-type sources located at 40 m ≤ y ≤ 60 m for t = 1000 years (M = number of terms summed for inverse generalized integral transform and N = number of terms summed for inverse finite Fourier cosine transform).When we investigate the required M for inverse generalized integral transform, N = 16 for the finite Fourier cosine transform inverse are used.When we investigate the required N for inverse finite Fourier cosine transform, M = 6400 for the generalized transform inverse are used.

Table 8 .
Transport parameters used for verification example 2 involving the ten-species transport problem used bySrinivasan and  Clement (2008b).

Table 10 .
Transport parameters used for example application involving the five-species dissolved chlorinated solvent problem used by BIOCHLOR.

Table 11 .
Coefficients of Bateman-type boundary source used for example application involving the five-species dissolved chlorinated solvent problem used by BIOCHLOR.