Lesson 6: Flow and Solute Transport Processes in 1D Systems
Lesson 6: Flow and Solute Transport Processes in 1D Systems hmg148Overview
This lesson introduces physical processes, including advection, diffusion, and dispersion processes, in 1D system. This is the first time in this course that we introduce the space dimension. An example will be shown about how to set up a one-dimentional (1D) flow and transport simulation in a homogeneous column in CrunchFlow.
Learning Outcomes
By the end of this lesson, you should be able to:
- Understand the definition of advective, diffusive, and dispersive transpor;
- Grasp the meaning of the one-dimentional advection-dispersion equation (ADE);
- Comprehend the important parameters and the application of dimensionless numbers (Péclet number) in unifying different observations;
- Understand what controls the shape of tracer breakthrough curves;
- Simulate 1D flow and transport of non-reactive tracers using CrunchFlow.
Lesson Roadmap
| Required Readings |
|
|---|---|
| Optional Reading (if phreeqc) |
|
| To Do |
|
Questions?
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
6.1 Background
6.1 Background mjg8Understanding flow and transport processes in the natural subsurface are important for a wide range of applications and disciplines. Flow and transport processes play a critical role in ground water and surface water management and environmental protection, energy extraction from deep hydrocarbon reservoirs, chemical weathering, and soil management. Here we primarily discuss fundamental flow and transport processes in natural subsurface systems.
For a conservative tracer that does not participate in reactions, advection, dispersion and diffusion are the major transport processes that control its transport. Advection, also called convection in some disciplines, determines how fast a tracer moves with the fluid flow; dispersion and diffusion processes are driven by concentration gradients and/or the spatial variation of flow velocities and therefore determine the extent of spreading.
6.2 Advection
6.2 Advection jls164Advection is the transport process where solutes flow with the bulk fluid phase. This is like you let go of yourself when you swim so you have the same velocity of the flowing fluid. The advective flux, of the solute, can be expressed as
(1)
where is the porosity of porous media; v is the linear fluid velocity in poroud media (m/s); and C is the solute concentration (mol/m3). Flow through a porous medium is described with Darcy’s Law:
(2)
where is the Darcy flux that is proportional to the gradient in the hydraulic head ; K is the hydraulic conductivity (m/s); One can also write Darcy’s Law in terms of hydraulic head by defining the hydraulic head as
(3)
where z is the depth (m), P is the fluid pressure (Pa), is the fluid density (kg/m3), and g is the acceleration of gravity (9.8 N/m2). The hydraulic conductivity (m/s) can be expressed as
(4)
where is the permeability of the porous media (m2) and is independent of fluid property, is the fluid hydraulic viscosity (Pa·s). Representative values of hydraulic conductivity and permeability are listed in Table 1 for various subsurface materials.
| K (m/s) | 100 | 10-1 | 10-2 | 10-3 | 10-4 | 10-5 | 10-6 | 10-7 | 10-8 | 10-9 | 10-10 | 10-11 | 10-12 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
(m2) | 10-7 | 10-8 | 10-9 | 10-10 | 10-11 | 10-12 | 10-13 | 10-14 | 10-15 | 10-16 | 10-17 | 10-18 | 10-19 |
| Unconsolidated Sand & Gravel | Clean Gravel | Clean Sand or Sand & Gravel | Very Fine Sand, Silt, Loess, Loam | ||||||||||
| Unconsolidated Clay & Organic | Peat | Stratified Clay | Unweathered Clay | ||||||||||
| Consolidated Rocks | Highly Fractured Rocks | Oil Rocks Reservoir | Sandstone | Limestone | Granite | ||||||||
By combing Eqn. (2)-(4), Darcy’s Law can also be written in terms of the fluid pressure, permeability, and the viscosity
(5)
Here, is the fluid pressure gradient. If the gravity term is negligible compared to the pressure gradient, Eqn. (5) can be simplified to
(6)
The characteristic time of the advection is the residence time, i.e., how long does a fluid particle stays within a given system. The residence time, , can be calculated as follows:
(7)
Here Vpore is the pore volume (m3), Q is the flow rate (m3/s), L is the length, A is the cross-section of the porous media in the direction perpendicular to the flow.
Darcy’s Law is applicable at the continuum scale where a representative elementary volume (REV) is significantly larger than the average grain size. The range for the validity of the Darcy’s Law can be checked using the Reynolds number Re:
(8)
where the volumetric flow rate Q (m3/s) and p is the perimeter of a channel or grain size (m). The upper limit of the validity of the Darcy’s Law is when Re is between 1 and 10 [Bear, 2013].
6.3 Diffusion
6.3 Diffusion jls164Molecular diffusion is driven by concentration gradient and is described by the Fick’s First Law:
Here Jdiff is the diffusive mass flux per unit area (mol/m2/s); D0 is the molecular diffusion coefficient in water (m2/s); and x is the distance (m).
Diffusion coefficients in water are typically in the order of 10-9 m2/s and depend on chemical species, temperature and fluid viscosity. Table 2 lists the diffusion coefficients in water at room temperature for some species. Given the diffusion coefficient at a specific temperature, the diffusion coefficient at another temperature can be calculated as follows:
where is the diffusion coefficient (m2/s) at a reference temperature T0 (K), and and are the dynamic viscosities at temperatures T and , T0 respectively.
| Cations | D0 (10-9 m2/s) | Anions | D0 (10-9 m2/s) |
|---|---|---|---|
| H+ | 9.311 | OH- | 5.273 |
| Li+ | 1.029 | F- | 1.475 |
| Na+ | 1.334 | Cl- | 2.032 |
| K+ | 1.957 | I- | 2.045 |
| NH4+ | 1.957 | NO3- | 1.902 |
| Mg2+ | 0.706 | HCO3- | 1.185 |
| Ca2+ | 0.792 | HSO4- | 1.331 |
| Al3+ | 0.541 | H2PO4- | 0.879 |
| Fe2+ | 0.719 | SO42- | 1.065 |
| Fe3+ | 0.604 | CO32- | 0.923 |
Fick’s second law combines the mass conservation and Fick’s first law:
where t is the time (s). Eqn. (9) can be analytically solved with appropriate boundary conditions.
Diffusion Coefficient De in Porous Media
Diffusion in porous media differs from diffusion in homogeneous aqueous solutions. Diffusion in porous media occurs through tortuous and irregularly shaped pores, as shown in Figure 1, and is therefore slower than that in homogeneous solutions. The effective diffusion coefficient describes diffusion in porous media and can be estimated using several forms:
Here F is the formation resistivity factor (dimensionless), an intrinsic property of porous media; the cementation factor (dimensionless) m quantifies the resistivity caused by the network of pores. Reported cementation factors vary between 1.3 and 5.0 [Brunet et al., 2013; Dullien, 1991]. For many subsurface materials, m is equal to 2 [Oelkers, 1996]. TL is the tortuosity (dimensionless) defined as the ratio of the path length in solution relative to the tortuous path length in porous media.

6.4 Dispersion
6.4 Dispersion jls164Dispersion describes the mixing of a solute due to fluctuations around the average velocity. This is caused by three factors: 1) microscopic heterogeneity, which make the fluid moves faster at the center of the pore and slower at the water grain boundary; 2) variations in pore sizes, in which cases fluid particles will move through larger pores faster; 3) variations in path length, causing some fluid particles going longer paths than others.
Mechanical Dispersion
Mechanical dispersion is a result of variations in flow velocities. Dispersion coefficients in porous media is typically defined as the product of the average fluid velocity and dispersivity :
where u is the average flow velocity (m/s) and refers to the dispersivity (m). In systems with more than one direction, the longitudinal dispersivity DL in the principle flow direction is typically higher than DT in the direction transverse to the main flow.
Dispersion is a scale-dependent process, with larger dispersivity values observed at larger spatial scales. At the column scale, a typical dispersivity may be on the order of centimeters. At the field scale, apparent dispersivities can vary from 10 to 100 m, as shown in Figure 2.
Hydrodynamic dispersion
The spreading of the solute mass as a result of diffusion and dispersion is similar to diffusion. This has led to the use of Fick’s First Law to describe the dispersion process as follows:
where Dh is the hydrodynamic dispersion coefficient defined as the sum of effective diffusion coefficient De and mechanical dispersion coefficient Dm:
As such, the hydrodynamic dispersion includes both diffusion and mechanical dispersion processes.
6.5 Advection Dispersion Equation (ADE)
6.5 Advection Dispersion Equation (ADE) jls164By combining the transport processes outlined above, we can derive an expression for the mass conservation of a non-reactive solute as follows:
(1)Substitution of Eqn. (1) and (13) into Eqn. (16) yields
(2)Equation (17) is the classical Advection-Dispersion equation (ADE). For one-dimensional systems, Equation (17) is simplified into
(3)Analytical solution of ADE is available for homogeneous porous media [Zheng and Bennett, 2002]:
(4)with the initial and boundary conditions:
(5)where C0 is the injecting concentration of the tracer, and erfc(B) is complementary error function:
(6) Here, erf(B) is error function, a special function of sigmoid shape that ocuurs in probability, statistics, and partial differential equations describing diffusion.6.6 Dimensionless numbers
6.6 Dimensionless numbers jls164Péclet number (Pe) is often used to describe the relative importance of advection and dispersion/diffusion in terms of their respective time scales :
(1) (2) (3)where L is the length of the domain of interest (m), u is the average Darcy flow velocity in the direction of interest (m/s), $D_{h}$ is the dynamic dispersion coefficient (m2/s). There are also some mathematical equations to define the time scales of these processes with similar concepts, mostly depending on the selected characteristic length [Elkhoury et al., 2013; Huysmans and Dassargues, 2005; Steefel and Maher, 2009; Szymczak and Ladd, 2009]. For example, L can also be the grid spacing (m) or correlation length (m) [Huysmans and Dassargues, 2005]. As shown in Figure 3, increasing Pe values indicate increasing dominance of advective transport and sharper front in breakthrough curves.
6.7 Setting Up 1D Flow and Transport Simulation in CrunchFlow
6.7 Setting Up 1D Flow and Transport Simulation in CrunchFlow jls164Please watch the following video: Advection-Dispersion Equation (ADE) for non-reactive tracers (16:42)
Advection-Dispersion Equation (ADE) for non-reactive tracers
So in is this video we'll be talking about advection-dispersion equation. This is first time actually we're not talking about chemistry but talking about the physical processes. So what do we have? Thinking about this is, let's say you have chemicals that are nonreactive. So the system you have we could use it as example is a column. So think about, for example, you are doing an experiment. You pack up a column, a column with like sand grains. And then you inject a chemical, let's say bromide, from the left into the right. So essentially you would have these chemical species you inject for a short period of time. You can imagine this is chemical will be moving along with the flow. So over time it will be eventually moving out of system with certain velocities. So v is the velocity here. And you are talking about here for example is the last of this column is L. So we want to know, as you can imagine, the constitution of this chemical would change with time and with space. Depends on what time you're looking at the snapshot. This chemical species might be here, might be here and other places there's no of this chemical species because we consider maybe, for example, we start with clean water.
So how do we mathematically solve this kind of equation and get the solution so we know the concentration of this chemical species in different time and different locations? And this is what we are going to talk about which is the advection-dispersion equation. We call it ADE with the other reaction terms. So ADE, when we think about it, I'm not going to derive in detail of where is this equation coming from. In general these equations are coming from these mass conservation principles. So if we look at these terms, different terms of the ADE, so the first home that we called the mass accumulation term. And it should have the units of accumulation, like mass per time. It's how fast things are changing. And this C is the concentration of the chemical species in the water phase. So C here is concentration of this-- let's call that tracer in water So everything we are solving is for how much they have in the water because that's what we really care about. And this should have the units of, for example, mass per volume.
So the first term is the mass accumulation. The second term here we call advective transport. So this is a process where there's a chemical almost you think about swimming. The chemical is a tracer. It actually flows together with the water at the same speed as the water flow. So that's called advective transport. And the last term is a dispersive transport. It's somewhat similar to diffusion process. It's driven by, first of all, concentration gradient, but also in this type of medias. You have grains. They tend to have different grain size, like different size of pore space and some spatial [INAUDIBLE] that leads to different flow velocity, mixing processes that lead to the concentration difference in different places. So these are three major terms in this equation. And for phi-- so this is phi is we call porosity, is how much space you have, how much pore space you have in a given volume. So the porous media has both pore space and solid space. So this is how much pore space in terms of percentage, per unit of total volume.
This v is velocity, flow velocity. It's a linear velocity in the pore space. So it's different from the u. We usually we call that's the velocity. And the relationship between the two, typically we knows that-- velocity-- we know that u will be equal to 5 porosity times the linear velocity. So there's a relationship between the two. And this Dh-- Dh is a very important term for this dispersive transport. This would be equal to-- I talk about it could be coming both from diffusion and the mechanical dispersion in porous media. So we have this kind of two terms adding together. This is coming from diffusion in porous media-- diffusion coefficient in porous media. And then the second term is essentially taking into account the mechanical dispersion. And you can think about alpha is called dispersivity, which is a parameter to describe how fast mechanical dispersion happens. And it's related to velocity. The whole term is also related to the velocity. So the flow, you actually have larger term of this.
So we have this equation. All the terms are here. And usually when this equation, there's an analytical solution for this equation if you give the right initial boundary condition. When we're numerically solving this, we will be discrete sizes of this equation in time and space, and then you get a solution for that. But before we do that, typical we need two conditions, right? This is the first time we introduced a space diminishing, x here. Before in the mineral dispersion precipitation, last time we actually induced the time. So you'll notice here that this equation, we call this is a partial differential equation, meaning it has two independent variables. One is the time. The other is space. So this is first time we introduce a space dimension here.
Now in order to solve this equation we need to know, like at t equal to 0, what are the concentrations? So this is initial condition. When t equal to 0, what are the concentrations? Now usually this is given for a given system. And we also need to know at the two boundaries, x equal to 0, x equal to L, what are the conditions that is specified? Is it, for example, a no-flow boundary or pure advective, or pure dispersive? These are different types of boundary conditions you can specify. But you will need to, in both the initial condition boundary conditions, to solve this. And different types of condition will give you different solutions because it matters what is the concentration in the t equal to 0. If your already start with something high you will see a different concentration versus, for example, at t equal to 0 you have clean water.
So terms of a solution, let's assume we have done all this work to solve the equation. What do we expect to see after we solve this equation? So I'm talking about a system. I'll be using an example system. So you will be injecting kind of shot pulse of this chemical bromide into a system. So if you conceptually think about it, at t equal to 0 it's somewhere here. And let's say everywhere it started with clean water. So you will have, at some point here, let's say you have kind of a little bit smeared at other time. And as time goes on you will have more and more smeared, but going in the direction of flow and moving along and maybe become more. So total mass will remain the same, but at the end you will see kind of wider and wider over time and over a longer distance. So this is conceptually how you would think about this solution you expect to see. Now when we think about from mathematical term, let's draw this. After we solve it, let's see. We'll look at the concentration as a function of distance. And what do you expect to see at different time?
So first of all, let's say at initial time you probably would see something like this. And here you see a pulse of this chemical. So is t at about 0, maybe a little bit past 0. But you think about, OK, over time, this pulse will move along. And then you should see different distributions. So this would be, let's say, at t equal to v times some small consolation, tu with C at another place. But this will become a little bit wider, smeared out. This is t1. And over time as you're going further distance, this becomes like more and more smooth. And I'm not drawing accurate. So total mass will remain the same. And I don't think I'm drawing nicely in terms of mass conservation. But in any case you will see over time total mass remain the same. But then the center of this will move along. So it's a v. So if it's a speed of this plume, how fast it moves over a certain time will be determined by-- so speed, v determines the speed of the center, of the rate of the center moving to downstream. But also another case will be-- we notice this Dh, which is a very important parameter as well. So if you're think about two different situations, let's say they both have the same v, but they will have different Dh values. Let's say we have another case with much higher Dh. Your will probably see something like this, a wider distribution. Still the same total mass, but it's much more spread. So this will be representing your large-- I'm sorry, a small Dh. The blue one would be representing a much larger Dh.
So now I draw the Dh. It will be more spread. So essentially you can think about Dh determines the width of the plume. So this is t1, t2. And over time, for example, if there's a t3 we care about, it should be even more spread out. This is t3. So if we wait long enough, this is going to spread a lot. So that's the type of solution you would expect to see when you have pulse of injection of a tracer. Now I'm just very briefly mentioning the characteristic time. There are several times we think is important. Once is residence time. This is directly related to how fast the flow goes and how long it actually stays in this column, how long this tracer actually stays in this column. So this we call tau a is equal to porosity times the length over mu. Or you can just say L over v. This is residence time. Another time is how long it takes for the dispersion process to uniform the whole concentration field. So this is what we call tau d. It's related to dispersion. So it's L squared divided by Dh. You can almost think about this as like if it's not open system, it would be how long it takes for diffusion to uniform, to homogenize the whole concentration field. And then a lot of times we define this Peclet number is tau d over tau a. So it's a relative between these two time scales or the rate of flow versus the rate of dispersion. And this should have a unit of L u phi Dh. This is a dimensionless number. So these times will determine like how fast-- you will realize this are kind of grouping dimensionless numbers. The Peclet number will determine the relative importance of dispersion versus advection, which in homework I asked you to do some exercises under different conditions with different Peclet numbers.
Example 6.1: Setting up 1D Flow and Transport Simulation in CrunchFlow
Click here for Example 6.1 CrunchFlow file package
This example introduces setting up numerical simulation of the flow and transport processes for a non-reactive tracer in a 1D column of 10 cm long. A tracer Br- is injected into the column at the concentration of 1.2×10-4 mol/L. The permeability of the column is 1.75×10-13 m2 and the porosity is 0.40. A constant differential pressure is imposed at the x direction and results in a Darcy flow velocity of 4.20×10-6 m/s. The molecular diffusion coefficient in water D0 for the tracer bromide is 1.8×10-9 m2/s (between that of F- and Cl- in Table 1). The cementation exponent m is 1.0. The dispersivity α is 0.07 cm. In order to keep consistent with the value of water viscosity provided in the CrunchTope original code, the water viscosity we applied here is 1.00×10-3 Pas at 20 0C.
- Before the numerical experiment:
- Calculate the pressure gradient for the Darcy flow velocity.
- Calculate the effective diffusion coefficient De, mechanical dispersion coefficient Dm, hydrodynamic dispersion coefficient Dh. Which term dominates in Dh?
- Calculate the characteristic time for advection (residence time) and for dispersion, and Pe number. Based on the Pe number, predict which term dominates during the transport process.
- Set up in CrunchFlow using 100 grid blocks and run the simulation:
- Plot the bromide breakthrough curve at the end of the column (C/C0, with C0 being the injection concentration) versus residence time.
Click for Solution
Click here for input and database file package
Solution 1:
For the 1-D flow, Darcy flow equation can be simplified to:
Thus,
This means that for a 10 cm (0.1 m) column, the pressure difference between one end to the other end of the column should be 2.40E3 Pa. Note that in Crunchflow, the pressure units is Pascal.
- Residence time:
Characteristic time for hydrodynamic dispersion:
CrunchFlow setup brief: This is the first time that we set up spatial dimension in CrunchFlow. Please read CrunchFlow manual the DISCRETIZATION key word block (page 46 – 47), INITIAL_CONDITION and BOUNDARY_CONDITION (page 72) as well as the TRANSPORT keyword block (page 72-75). The CrunchFlow exercise example 6 also teaches how to set up non-reactive tracer test.
Solution 2:

Figure 4. Breakthrough outlet concentration of Br- as a function of number of residence time. C0 is the initial injecting concentration of the tracer Br- while C is the tracer concentration at the outlet.
6.8 Homework Assignment
6.8 Homework Assignment jls164Assignment
- (60 points) As an extension of the in-class example, we will do another tracer test that only inject the tracer for the first 1/10 of the residence time. The rest of the time we will inject clean water into the column. Please note that in order to do this you will need to generate two input files with one inputfile having the keywords "restart", "save_restart", and "later_inputfile". Please read the relevant section in CrunchFlow manual for these keywords. The exercise 12 in the exercise folder has sample input files on using these key words.
- 1.1 Please plot the spatial profile of the tracer (concentration as a function of distance) at 1/10, 2/10, 4/10, and 8/10 of the residence time. Please plot all spatial profiles in one figure. Describe how the shapes of the Br concentration curve change over time. What causes that change?
- 1.2 Change the dispersvity to 0.1 and 10.0 of the base value in example 6.1. Rerun these two additional simulations. And do the same plot as in question 1). Observe the difference in the shape of Br- concentration spatial profile with different values. What are the effects of changing values?
- (optional) (60 bonus points) As an extension of the in class example with continuous injection of Br- from the inlet, we will do a sensitivity analysis to understand the effect of parameters on the breakthrough curves in the same 1D column, including flow velocity, dispersivity, molecular diffusion coefficient, and cementation factor. In each sub question below, only one parameter or condition is changed. For each condition, please
- 2.1 calculate the characteristic times (residence time and dispersion time) and Pe;
- 2.2 draw the Br breakthrough curves (as a function of residence time) by changing only the parameter discussed in each sub-question, with all other parameters remaining the same as those in example 1. Note that each middle number is the number used in example 6.1 so you do not need to run the simulation again.
- flow velocity: 0.5, 1.0, and 2.0 of the flow velocity in example 6.1. Plot all three BTC in one figure and describe the effect of flow velocity on breakthrough curves (BTCs).
- dispersvity : 0.1, 1.0, and 10.0 of the base value in example 6.1. Plot all three BTCs in one figure and describe the effects of dispersivity on BTC.
- Molecular diffusion D0: 0.1, 1.0, and 10.0 of the base value in example 6.1. Plot all three BTCs in one figure and describe the effect of D0 on BTC.
- cementation factor m: 1.0, 3.0, 5.0 of the base value in example 6.1. Plot all three BTC in one figure and describe the effect of cementation factor on BTC.
- compile all previous cases in 1) – 4) and make a table that lists values of residence time , dispersion time , De, Dm, Dh, v, and Pe. Plot all breakthrough curves together, label the Pe values for different cases in different ranges. What do you observe and why?
6.9 Summary and Final Tasks
6.9 Summary and Final Tasks hmg148Summary
In this lesson, we covered definition and principles of transport processes including advection, diffusion, and dispersion processes. We also discussed the Advection-Dispersion-Equation (ADE), its analytical solution, and how to set up solving ADE in CrunchFlow.
Other Relevant Textbooks:
-
Flow and Transport in Porous Formations (1989) by Gedeon Dagan;
-
Dynamics of fluids in porous media (2013) by Jacob Bear, Courier Dover Publications;
-
Applied contaminant transport modeling (2002), 2nd edition, by Chunmiao Zheng and Gordon D. Bennett, John Wiley and Sons, Inc.
References
References hmg148Bear, J. (2013), Dynamics of fluids in porous media, Courier Dover Publications.
Brunet, J. P. L., L. Li, Z. T. Karpyn, B. G. Kutchko, B. Strazisar, and G. Bromhal (2013), Dynamic Evolution of Cement Composition and Transport Properties under Conditions Relevant to Geological Carbon Sequestration, Energy & Fuels, 27(8), 4208-4220.
Dullien, F. A. (1991), Porous media: fluid transport and pore structure, Academic press.
Elkhoury, J. E., P. Ameli, and R. L. Detwiler (2013), Dissolution and deformation in fractured carbonates caused by flow of CO< sub> 2</sub>-rich brine under reservoir conditions, International Journal of Greenhouse Gas Control, 16, S203-S215.
Gelhar, L. W. (1986), Stochastic subsurface hydrology from theory to applications, Water Resour Res, 22(9S), 135S-145S.
Haynes, W. M. (2012), CRC handbook of chemistry and physics, CRC press.
Huysmans, M., and A. Dassargues (2005), Review of the use of Péclet numbers to determine the relative importance of advection and diffusion in low permeability environments, Hydrogeol J, 13(5-6), 895-904.
Oelkers, E. H. (1996), Physical and chemical properties of rocks and fluids for chemical mass transport calculations, Reviews in Mineralogy and Geochemistry, 34(1), 131-191.
Steefel, C. I., and K. Maher (2009), Fluid-Rock Interaction: A Reactive Transport Approach, Rev Mineral Geochem, 70, 485-532.
Steefel, C. I., D. J. DePaolo, and P. C. Lichtner (2005), Reactive transport modeling: An essential tool and a new research approach for the Earth sciences, Earth Planet Sc Lett, 240(3-4), 539-558.
Szymczak, P., and A. J. C. Ladd (2009), Wormhole formation in dissolving fractures, J Geophys Res-Sol Ea, 114.
Zheng, C., and G. D. Bennett (2002), Applied contaminant transport modeling: theory and practice, 2nd ed., 621 pp., John Wiley and Sons, Inc., New York.