# 2.2: Geological Applications of the Diffusion Equation

- Page ID
- 3527

## General Diffusion Equation

**Diffusivity: **Note that \(\kappa\) has units of \(\frac{m^2}{s}\), which does not depend on temperature or heat flow. Diffusivity, in general, describes the spreading out of some quantity (e.g. temperature) over space in some unit of time. Other diffusion processes will have other diffusivities, but they all will have the same unit. For each of the equations below, the solution is found by defining an initial condition (the initial temperature, topography, or concentration at every position) and two boundary conditions (the value of the temperature, topography, or concentration at two positions...they can be fixed or functions of time). We will go over these conditions in more detail in the next few sections.

**1-D vs 2-D/3-D**: Heat spreads out in 2-D or 3-D because of the addition of the \(\frac{\partial^2 T}{\partial y^2}+\frac{\partial^2 T}{\partial z^2}\) term to the diffusion equation

\[\frac{dT}{dt}=\kappa\frac{d^2T}{dz^2}.\]

If \(\kappa\) is a constant, it will work the same in 1-D as it does in 2-D or 3-D, heat will just spread out in multiple directions.

The contours on the above figure are analogous to the contours on a topographic map. The contours are spaced closely in the x direction and so the temperature gradient is large there. Thus, the higher the degree of curvature on the graph indicates the direction of maximum cooling/heating.

**Cooling of the Oceanic Lithosphere**

\[\frac{dT}{dt}=\kappa\frac{d^2T}{dz^2}\]

where \(\kappa\) is the thermal diffusivity. The solution for this equation is:

\[T(x, t)=(T_s-T_m)erfc(\frac{y}{2\sqrt{\kappa t}})+T_m\]

where \(T_s\) is the surface temperature, \(T_m\) is the mantle temperature, and y is the depth into the lithosphere. We will use this solution all the time, as it indicates the thermal structure of the oceanic lithosphere. It does not take into account reheating due to volcanism or mechanical feedback, but it is used to define the initial condition in models. We will use it to calculate the ridge push force.

As you can see above, in order to solve the partial differential diffusion equation, we took several large mathematical steps. Let's explore how we got there in more detail. Solving partial differential equations (PDE) might seem challenging at first, especially if you have never taken a class on PDE before, but by breaking it into a few simple steps, it becomes less daunting. Let's start with the basic diffusion equation:

\[\frac{dT}{dt}=\kappa\frac{d^2T}{dz^2}\]

Now we will integrate to find \(T=f(x, t)\). \(\frac{dT}{dt}\) on takes the derivative in terms of t and holds x constant. For __each__ derivative there will be a constant of integration.

\[T(x)=\int_{0}^{x_o}x dx=\frac{1}{2}x^2+C_o\]

First, In order to find \(C_o\), we need a boundary condition. How many boundary conditions do we need? One for every derivative. If we have \(\frac{d}{dt}\), we need one initial condition. \(\int dt\rightarrow C_o \rightarrow time\)⇒T at all x locations at t=0. In more detail, for the 1-D diffusion equation this means:

- One condition on time: this is called the initial condition. Simply put, what is the starting value for the quantity you are calculating at every point along your profile? This might be a constant or some function of position.
- Two conditions on space: these are called boundary conditions. What is the value for the quantity you are calculating at each end of your profile (or at infinity depending on the specific problem) for ALL times? These values might be constants or some function of time.
- If we have \(\frac{d^2}{dx^2}\), we need two boundary conditions. \(\int(\int dx)dx \rightarrow C_1x+C_2\rightarrow\) space and distance ⇒T at all times at x=0.

Second, we need to know values of the physical quantities that show up in our equation. For the diffusion equation, the only constant is the diffusivity that is relevant for the particular physics.

Third, we need to either know how to solve a PDE or where to look up the general solution. For almost every diffusion equation you come upon, the solution to the equation will include dependence on the error function.

\[erf(\eta)=\frac{2}{\sqrt{\pi}} \int_{0}^{\eta} e^{-\eta^{\prime 2}} d\eta\prime\]

The parameter \(\eta\) depends on the distance, time, and the diffusion coefficient, \(\eta=\frac{z}{2\sqrt{\kappa t}}\). Note that if you set your \(\eta=1\), you are back to the simple equation that we used to determine the length-scale and time-scale for diffusion. The value of the error function for any particular value of \(\eta\) can be calculated using a built-in function or looked up in a table. There is also a related function called the complementary error function.

\[erfc(\eta)=1-erf(\eta)\]

Here is a plot of the error function and the complementary error function (in the plot x is \(\eta\)). The final step to solve the PDE is to use the initial conditions and boundary conditions to determine the constants of integration in the general solution.

Now that we have a general idea of how to solve a PDE, let's return to the oceanic lithosphere example.

Example \(\PageIndex{1}\): Cooling of a Lithospheric Plate

Lithosphere forms by cooling of the mantle as the tectonic plate moves away from a spreading ridge. At the ridge, the mantle temperature of \(T_m\) almost reaches the surface (z=0). As the material formed at the ridge moves away from the spreading center at a rate equal to the speed of the tectonic plate (\(u\)), the mantle cools from above to become the lithosphere. The temperature at the surface of the lithosphere, \(T_s\), is essentially constant because the surface is at the bottom of the ocean (\(T=0^oC)\). The age of the plate at a distance \(x\) from the spreading center is determined by the speed of the plate: \(t=\frac{x}{u}\). A typical value for the thermal diffusion coefficient (thermal diffusivity) is \(\kappa_T=1.0x10^{-6}\frac{m^2}{s}\).

*What is the initial condition for the temperature?*

At time t=0, the temperature is equal to the mantle temperature at all positions in z, \(T(z, t=0)=T_m\). This is even true at the surface. The temperature at the surface changes to zero at the instant of time just after t=0.

*What are the boundary conditions for the temperature?*

The temperature far away from the surface will stay at the mantle temperature, T_{m}, and the surface temperature will stay at T_{s}=0. At z=0, T(0,t)=T_{s}. As \(z\rightarrow\infty\), T=T_{m}. Put in another way:

at t=0 T(z=0)=\(T_m\) and, at z=0 T(t>0)=\(T_s\)

t=\(0_+\) T(z=0)=\(T_s\) z=\(\infty\) T(t>0)=\(T_m\)

*What is \(\eta\) in the error function equation?*

\[\eta=\frac{x}{2\sqrt{\kappa t}}\]

*Solution*:

\[T(x, t)=(T_s-T_m) erfc(\frac{x}{2\sqrt{\kappa t}})+T_m\]

Let's try another example.

Example \(\PageIndex{2}\): Fault Scarp Erosion

The change in shape of a fault scarp profile as a function of time due to erosional processes is often modeled using the diffusion equation. In this case the topography (mass) **diffuses** (erodes) from the uplifted block to the down-dropped block. For the case of an initial vertical offset, as the diffusive erosion progresses the sharp discontinuity in surface topography will become smoother with the slope across the fault with decreasing time.

For this example, we will consider a 1-D topography profile across a fault with total vertical offset \(\Delta H(x=0)=2a\), and \(H_{x_{max}}=a\) and \(H_{x_{min}}=-a\). The fault could occur in a region with an overall background slope given by b (units of \(\frac{m}{km}\); b=0 is a flat surface). The diffusion coefficient for topography will depend on the type of rock forming the fault and surface processes. Typical diffusion coefficient values are: \(\kappa_H=0.5-1.0 \frac{m^2}{kyr}\), which is \(\kappa_H=1.5-3.0\) x10^{-11} \(\frac{m^2}{s}\).

*What is the initial condition for topography?*

At time t=0, the fault has just formed so the topography on either side of the fault has a specified value. For a flat background surface:

At t=0 and \(x\leq 0, H(x)=a\) or \(H(x)=a+bx\)

At t=0 and \(x>0, H(x)=-a\) or \(H(x)=-a-bx\)

*What are the boundary conditions for the topography?*

For the fault scarp diffusion problem, erosion will only modify the topography in the region near the fault. Far away from the fault, the total offset will not change from the initial value:

As \(x\rightarrow -\infty\), for x\(\leq 0\), \(H(x, t)\rightarrow a+bx\)

As \(x\rightarrow \infty\), for x\(> 0\), \(H(x, t)\rightarrow -a-bx\)

*What is \(\eta\) in the error function equation?*

\[\eta=\frac{x}{2\sqrt{\kappa_H t}}\]

For the case of erosion of a fault scarp with initial height, H, the diffusion equation can be written as:

\[\frac{\partial H}{\partial t}=\kappa_H(\frac{\partial^2H}{\partial x^2})\]

where \(\kappa_H\) is the diffusivity for fault scarp topography.

*Solution*:

Using the initial and boundary conditions the diffusion equation can be solved for the change in shape of the topography as a function of time and position:

\[H(x, t)=erf(\frac{x}{2\sqrt{\kappa_H t}})\alpha+bx\]

where \(\alpha\) is half the total fault offset and \(b\) is the background slope of topography.

Let's do one last example to practice our new PDE solving skills.

Example \(\PageIndex{3}\): Diffusion of Chemical Species (a Solute in Water)

Contamination of groundwater by toxic spills can be modeled as a diffusive process, where the concentration of a toxic substance diffuses out away from the source due to a chemical gradient. Basically, the toxin wants to move form a region of high concentration to a region of low concentration in order to minimize the total energy of the system (thermodynamics concept).

For this example imagine there is a gasoline station with a leaking subsurface gasoline storage compartment. The rate of leakage is constant, so at the position of the leak there is a constant input value of toxin per volume, or constant concentration, C_{s} (s-spill). Far away from the leak, the concentration is zero (or in general some background value, C_{o}, depending on the solute of interest). For a 1-D profile, with the leakage located at x=0, the diffusion equation can be used to determine how the concentration of the toxin changes with time and distance away from the leakage. The diffusion coefficient for toxins will depend on the type of toxin, the medium it is diffusing through (groundwater in this example) and temperature. Typical concentration diffusion coefficients (\(\kappa_C\)), in groundwater at \(T=25^oC\), are: MTBE (Methyl tert-butyl ether, in gasoline) 9.4x10^{-10}\(\frac{m^2}{s}\), Benzene (in gasoline) 1.1x10^{-9}\(\frac{m^2}{s}\), and Methanol 1.65x10^{-9} \(\frac{m^2}{s}\).

*What is the initial condition for the solute concentration?*

At time t=0, the leakage starts, so the concentration at x=0 is C=C_{s}. Elsewhere the concentration is equal to the background value. But also note that at all times C=C_{s} since the spill continues at a constant rate.

*What are the boundary conditions for the solute concentration?*

Far away from the leakage, the concentration will be unaffected by the leakage: As \(x\rightarrow\pm\infty\), C=C_{o} (this is in fact two boundary conditions at plus and minus infinity). However, in addition, at the location of the spill the concentration stays constant because the leakage continues at a constant rate: C=C_{s} at x=0 for all times.

*What is \(\eta\) in the error function equation?*

\[\eta=\frac{x}{2\sqrt{\kappa_C t}}\]

For the case of diffusion of some chemical element denoted by C, the diffusion equation can be written as:

\[\frac{\partial C}{\partial t}=\kappa_C(\frac{\partial^2C}{\partial x^2})\]

where \(\kappa_C\) is the diffusivity of the chemical element and \(C\) is the concentration of the chemical element.

*Solution*:

Using the initial and boundary conditions the diffusion equation can be solved for the change in concentration of the toxin in the groundwater as a function of time and position

\[C(x, t)=(C_s-C_o)erfc(\frac{x}{2\sqrt{\kappa_C t}})+C_o\]

Error Functions

A note on error function: It is oftentimes simpler to use the error function instead of the complementary error function, as the error function starts at 0 and increases, like temperature.

Here is how we convert from erfc to erf:

For all times at z=0 \(T=T_s\) and, at t=0 T(z)=\(T_m\)

z=\(\infty\) \(T=T_m\)

\[T(z, t)=(T_s-T_m) erfc(\frac{z}{2\sqrt{\kappa t}})+T_m\]

at z=0, \(\eta\)=0, erfc=1 ⇒T_{s}-~~T~~+_{m}~~T~~⇒T=T_{o}_{s}

at z=\(\infty\), \(\eta=\infty\), erfc=0⇒T=T_{m}

at t=0, \(\eta=\infty\), erfc=0, ⇒T=T_{m}

We can rewrite this using erf(\(\eta\)).

\[T=(T_s-T_m)erfc(\eta)+T_m\]

Plug in

\[ \begin{align*} erf(\eta) &=(1-erf(\eta)) \\[4pt] &=(T_s-T_m)(1-erf(\eta))+T_m) \\[4pt] &=(T_s-T_m)-(T_s-T_m)erf(\eta)+T_m \\[4pt] &=T_s-(T_s-T_m)erf(\eta) \\[4pt] &=T_s+(T_m-T_s)erf(\eta) \end{align*}\]

\[T(z, t)=T_s+(T_m-T_s)erf(\frac{z}{2\sqrt{\kappa t}})\]