Lesson 3 - Understanding the Diffusion Equation

Before we discuss how we can model diffusion in GoldSim, it is necessary to first discuss how GoldSim represents diffusion mathematically. As we have mentioned several times, as a general rule, you should not use “black boxes” to build models. If you don’t understand how a tool works, you should not use it. Hence, in this Lesson we again will briefly discuss the actual equations that are used in GoldSim, and will also discuss (in general terms) how GoldSim solves the equations numerically. An understanding of these topics is necessary in order to properly model diffusive transport in GoldSim.

Let’s consider a tube example similar to what we discussed in the previous Lesson. This one is filled with sand and is full of water (i.e., the sand is saturated). Let’s assume we inject a slug of mass (uniformly over the entire cross-sectional area) at one end of a tube, but do not allow any water to flow through the tube at all. We will assume that the end of the tube where we inject the mass is closed (no mass can leave through that end).  However, at the other end of the tube there is a screen (that holds the sand in the tube) and the tube is connected to a small tank filled with water that is rapidly flushed (effectively maintaining the concentration in the tank at zero). What would happen to the concentration of the chemical in the tube?

What would happen, of course, is that the molecules would diffuse through the tube (due to the random thermal motion – Brownian motion - of the molecules) until eventually they had all diffused out into the tank (which is constantly flushed).  Diffusion is the macroscopic manifestation of the Brownian motion at the microscopic level. These random movements can be described statistically, and the evolution of the concentration in the tube can be described using Fick’s Second Law:

D_m is the molecular diffusion coefficient (dimensions of L2/T).

Due to the form of this equation, the nature of diffusion is that the diffusing substance tends to spread steadily from regions of high concentration to regions of lower concentration.

Molecular diffusion through water is a slow process (diffusion of the length of a meter through water would take on the order of 10 years for a typical molecule; while diffusion through air over the same distance would take on the order of a day). Nevertheless, this process can be critical for some applications, particularly those where engineered diffusive barriers are part of the design and the models have time scales of interest that are very large (e.g., radioactive waste management applications).

The molecular diffusion coefficient depends on the medium through which the chemical is diffusing (e.g., water or air), the temperature, and the size of the diffusing molecules.  In many mass transport models, however, we will find that the diffusion is taking place through a porous medium (such as in our example, in which the chemical is diffusing through a saturated sand). The diffusion occurs not in the sand particles, of course, but rather in the interstitial water that exists between the sand particles. In this case, the equation above needs to be adjusted for the properties of the porous medium:

In this equation, n is the porosity of the solid, and τ is the tortuosity. The tortuosity quantifies the degree to which diffusion is slowed down due to the tortuous (longer) path that the diffusing chemical must travel as it moves through the pore space.  We will discuss this input more in the next Lesson.

In GoldSim, we will simulate diffusion by using Cells to spatially discretize the system through which diffusion is occurring.  Each Cell will be well-mixed, and we will model a diffusive mass flux between the Cells based on the concentrations in each Cell, geometric properties, the diffusion coefficient, and the properties of the porous medium (if present) through which the chemical is diffusing.

To represent the system in this way, we need to approximate both the time derivative and the spatial derivative using finite differences. The (implicit, centrally-weighted) finite difference approximation for Fick’s Second Law for diffusion along a tube looks like this:

In this equation, we have assumed that there is no porous medium (if a porous medium was present, D_m would be multiplied by τ and n). In this equation, the tube has been discretized spatially into a fixed number of finite volumes of length Δx (where Δx is equal to L/N, with L being the total length of the tube and N is the number of Cells used to represent it). The equation involves the concentration in three adjacent volumes: j, j-1 (which is immediately to the left of j) and j+1 (which is immediately to the right of j). The superscript k refers to the timestep (k is the current timestep k-1 is the previous timestep, and Δt is the timestep length). There is an equation for every volume so, for example, if the tube was divided into 10 finite volumes, there would be a system of 10 equations with 10 unknowns for each contaminant being modeled. 

As we noted in the previous Lesson, discretizing the system in this way will produce some numerical errors.  Hence, as we will discuss later in this Unit, we will need to select our discretization in such a way that these errors are minimized. 

For now, however, let’s look at this equation a little closer so we can understand the implications of its general form. You will recall that when discussing advective transport between Cells, we wrote our governing equation in terms of mass rather than concentration (as this is how the equations are actually represented in GoldSim). In particular, if we had three Cells in series, the approximation for the rate of change of mass in the “middle” Cell (j) due to advective transport from an upstream Cell j-1 (and to a downstream Cell j+1) could be written as follows:

The first term on the right-hand side of the equation represents the mass transfer rate (into j) from the upstream Cell and the second term represents the mass transfer rate (out of j) to the downstream Cell.

It would be nice to write our diffusion equation in a similar way so we can compare it to how we represented advection.  To do this, we simply need to recognize that assuming no porous medium is present, C = M/V, where M is the mass in the Cell and V is the volume. We also know that V = A Δx.  Substituting these, the finite difference approximation for Fick’s Second Law for diffusion  looks like this:

Finally, let’s rearrange the right-hand side of the equation into two separate terms:

The first term on the right-hand side of the equation represents the mass transfer rate (into j) from the Cell to the left and the second term represents the mass transfer rate (into j) from the Cell to the right. In GoldSim, we refer to the (D_m A / Δx) term as the diffusive conductance. As noted previously, if we were simulating diffusion through a porous medium, the diffusive conductance would also include terms for the tortuosity and the porosity: (τ n D_m A / Δx).

Comparing this equation to the finite difference approximation for pure advection between Cells, we will note that the Q term for advection is analogous to the D_mA/ Δx term for diffusion (in fact, they have the same dimensions). You will recall that to define advection between Cells, we created an advective mass flux link (using the Outflows tab of the Cell element). When we did so, we specified the flow rate (Q).  In the next Lesson, we will discuss how we create a diffusive mass flux link between Cells in GoldSim. As you can guess, when doing so, we are going to need to specify τ, n, D_m, A and Δx.

One other point becomes clear when we look at how the advective mass transfer rate is represented compared to how the diffusive mass transfer rate is represented. Whereas the mass transfer rate from one Cell to another for advection depends only on the concentration of the upstream Cell, the mass transfer rate from one Cell to another for diffusion depends on the concentrations in both Cells (their difference). In fact, the mass transfer rate can reverse during a simulation depending on which of the two concentrations are greater. That is, during a simulation, the net mass transfer rate could potentially be positive for some parts of the simulation (e.g., mass moves from j+1 to j) and negative for other parts of the simulation (mass moves from j to j+1).