Ch12.4 Lake Pollution Revisited

• There are several processes by which pollution travels in lake, such as advection, dispersion, diffusion, and reaction.
• In Chapter 2.5, we assumed the lake was well mixed and that pollution was transported by water flow.
• This type of transport is advection, and is equivalent to convection of heat.

• We will assume pollution levels are different throughout lake, and that advection is responsible for movement of pollution.
• A PDE for pollution concentration will be derived, and solved analytically using the method of characteristics.
• Numerical solutions can be obtained using Math 366 methods, and software packages such as Maple.

Imagine if Americans switched from pounds to kilograms overnight.

There would be mass confusion!

\[ \small{ \frac{\partial C}{\partial t} = - v \frac{\partial C}{\partial x} + D_x \frac{\partial^2 C}{\partial x^2} + D_y \frac{\partial^2 C}{\partial y^2} + D_z \frac{\partial^2 C}{\partial z^2} - \lambda RC } \]

• Left side: rate of concentration change per unit volume.
• Right side: advective flux in x direction, solute flux due to dispersion in \( x, y, z \) directions, and biotic & abiotic decay.
• www.enviro.wiki/groundwaterflow

\[ \small{ \frac{\partial C}{\partial t} = - v \frac{\partial C}{\partial x} + D_x \frac{\partial^2 C}{\partial x^2} + D_y \frac{\partial^2 C}{\partial y^2} + D_z \frac{\partial^2 C}{\partial z^2} - \lambda RC } \]

Parameters

• \( v \) = advective transport velocity in \( x \) direction \( (L/T) \)
• \( D_x, D_y, D_z \) = hydrodynamic dispersion coefficients in \( x, y, z \) directions \( (L^2/T) \)
• \( \lambda \) = effective first order decay rate due to combined biotic and abiotic processes \( (1/T) \)
• \( R \) = equilibrium-based (reaction) reduction factor

\[ \small{ \frac{\partial C}{\partial t} = - v \frac{\partial C}{\partial x} + D_x \frac{\partial^2 C}{\partial x^2} + D_y \frac{\partial^2 C}{\partial y^2} + D_z \frac{\partial^2 C}{\partial z^2} - \lambda RC } \]

\[ \small{ \frac{\partial C}{\partial t} = - v \frac{\partial C}{\partial x} + D_x \frac{\partial^2 C}{\partial x^2} + D_y \frac{\partial^2 C}{\partial y^2} + D_z \frac{\partial^2 C}{\partial z^2} - \lambda RC } \]

In many lakes there is a channel of flow from entrance to exit that is rectangular in shape, where advection is main way that pollution travels. Thus we start by assuming the following:

• Pollution source will be from a feeder river at one end.
• Flow of water into and out of lake are the same.
• Lake is rectangular with constant volume \( V \).

• Pollution travels with flow of water (advection) only.
• Unidirectional flow of water from entrance to exit of lake.
• If \( x \) is in direction of water flow, and \( y \) and \( z \) are the other axes in 3D plane, then there is no change in pollution concentration in these other directions.
• Pollution is not in uniform concentration throughout lake.

\[ \small{ \begin{aligned} \frac{\partial C}{\partial t} &= - v \frac{\partial C}{\partial x} + D_x \frac{\partial^2 C}{\partial x^2} + D_y \frac{\partial^2 C}{\partial y^2} + D_z \frac{\partial^2 C}{\partial z^2} - \lambda RC \\ \frac{\partial C}{\partial t} & = - v \frac{\partial C}{\partial x} \end{aligned} } \]

• Let \( x \) be horizontal distance of pollution from entry point.
• Let \( C \) denote concentration of pollution at \( x \).
• The concentration \( C(x,t) \) of pollutant depends on both position \( x \) in lake and time \( t \).

• Let \( J(x, t) \) = pollution mass flux per unit area \( A \) per unit time at position \( x \) and time \( t \).
• Let \( F(x, t) \) = flow rate (\( V \)/time) across area \( A \).
  • \( F(a, t) \) = flow rate across area \( A \) at entrance
  • \( F(b, t) \) = flow rate across area \( A \) leaving lake.

Applied to total mass of pollutant in section of width \( \Delta x \),

\[ \small{ \begin{Bmatrix}\mathrm{rate \, of \, change\, of} \\ \mathrm{pollutant \, mass} \\ \mathrm{in \, section}\end{Bmatrix} = \begin{Bmatrix} \mathrm{rate \, of \, pollutant} \\ \mathrm{mass \, entering} \\ \mathrm{section \, at \,} x \end{Bmatrix} -\begin{Bmatrix} \mathrm{rate \, of \, pollutant} \\ \mathrm{mass \, leaving \, section } \\ \mathrm{at \, }x+\Delta x\end{Bmatrix} } \]

\[ \small{ \begin{Bmatrix}\mathrm{rate \, of \, change\, of} \\ \mathrm{pollutant \, mass} \\ \mathrm{in \, section}\end{Bmatrix} = \begin{Bmatrix} \mathrm{rate \, of \, pollutant} \\ \mathrm{mass \, entering} \\ \mathrm{section \, at \,} x \end{Bmatrix} -\begin{Bmatrix} \mathrm{rate \, of \, pollutant} \\ \mathrm{mass \, leaving \, section } \\ \mathrm{at \, }x+\Delta x\end{Bmatrix} } \]

\[ \small{ \begin{aligned} A \Delta x \frac{ \partial C }{ \partial t } &= J(x, t)A - J(x + \Delta x, t)A \\ \\ \lim_{\Delta x \rightarrow 0} \frac{ \partial C }{ \partial t } &= - \lim_{\Delta x \rightarrow 0} \frac{ J(x + \Delta x, t) - J(x, t) }{ \Delta x } \\ \\ \frac{ \partial C }{ \partial t } &= - \frac{ \partial J }{ \partial x } \end{aligned}} \]

From the previous slide,

\[ \small{ \frac{ \partial C }{ \partial t } = - \frac{ \partial J }{ \partial x } = -\frac{F}{A}\frac{\partial C}{\partial x}} \]

In above equation, we used \( J = (F/A)C \), where

\[ \small{ \left[J\right] = \frac{kg}{A \cdot sec} , \,\, \left[\frac{F}{A}C\right] = \left(\frac{\frac{V}{sec}}{A}\right) \frac{kg}{V} = \frac{kg}{A \cdot sec}} \]

Thus \[ \small{ \frac{\partial C}{\partial t} = - v \frac{\partial C}{\partial x}, \,\,\, v = \frac{F}{A},\,\,\, [v] = \left[\frac{\frac{V}{sec}}{A }\right] = \frac{L_x}{sec} } \]

Boundary condition at \( x=a \), where pollution enters lake:

\[ \small{ C(a,t) = \frac{g(t)}{F} } \]

Unit check:

\[ \small{ \left[\frac{g(t)}{F}\right] = \frac{\frac{kg}{sec}}{\frac{V}{sec}} = \frac{kg}{V} } \]

We also assume an initial pollution concentration each point \( x \):

\[ \small{ C(x,0) = P(x) } \]

Our first order PDE boundary and initial value problem:

\[ \small{ \begin{aligned} \frac{\partial C}{\partial t} &= - v \frac{\partial C}{\partial x} \\ C(a,t) & = \frac{g(t)}{F} \\ C(x,0) & = P(x) \end{aligned}} \]

Figure: Blue curve = \( P(x) \), orange curve = \( g(t) \).

• Pollution concentration is constant on characteristic lines

\[ \small{x = vt + K } \]

• The value of \( C(x,t) \) on boundaries will therefore propagate unchanged along these lines.
• Use this result to build analytical solution for \( C(x,t) \).

Thus for \( (x,t) \) on this characteristic line,

\[ \small{ \begin{aligned} t &= \frac{A}{F} (x - a) + K_1 \\ C(x,t) & = \frac{g(K_1)}{F} = \frac{g \left(t - \frac{A}{F} (x - a) \right)}{F} \end{aligned}} \]

Concentration of pollutant for \( \, a \leq x \leq b, \, t \geq 0 \, \):

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} P\left(x-\frac{F}{A}t\right), & \text{if } t<\frac{A}{F}(x-a) \\ g\left(t-\frac{A}{F}(x-a)\right), & \text{if } t>\frac{A}{F}(x-a) \end{array} \right. \\ } \]

Suppose initial and boundary conditions are given by

\[ \small{ \begin{aligned} C(x,0) & = P(x) = \frac{1}{4}(x-2)^2, \,\, 0 \leq x \leq 4\\ C(0,t) &= g(t) = e^{-(4-t)^2}, \,\, t \geq 0 \end{aligned}} \]

With \( v = 1 \), \( a = 0, b= 4 \), the solution

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} P\left(x-vt\right), & \text{if } t<\frac{x-a}{v} \\ g\left(t-\frac{x-a}{v}\right), & \text{if } t>\frac{x-a}{v} \end{array} \right.} \]

becomes

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} \frac{1}{4}\left(x-t-2\right)^2, & \text{if } t< x \\ e^{-(4-(t-x)^2)}, & \text{if } t > x \end{array} \right.} \]

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} \frac{1}{4}\left(x-t-2\right)^2, & \text{if } t< x \\ e^{-(4-(t-x)^2)}, & \text{if } t > x \end{array} \right.} \]

• CalcPlot3D Example

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} \frac{1}{4}\left(x-t-2\right)^2, & \text{if } t< x \\ e^{-(4-(t-x)^2)}, & \text{if } t > x \end{array} \right.} \]

• CalcPlot3D Example

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} \frac{1}{4}\left(x-t-2\right)^2, & \text{if } t< x \\ e^{-(4-(t-x)^2)}, & \text{if } t > x \end{array} \right.} \]

• Concentration as level curves over characteristics.

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} \frac{1}{4}\left(x-t-2\right)^2, & \text{if } t< x \\ e^{-(4-(t-x)^2)}, & \text{if } t > x \end{array} \right.} \]

Suppose only freshwater enters lake, with \( g(t) = 0 \).

\[ \small{ C(x,t)=\left\{ \begin{array}{ll} 0, & \text{if } t>\frac{x-a}{v} \\ P\left(x-vt\right), & \text{if } t<\frac{x-a}{v} \end{array} \right.} \]

For \( v=1 \), lake is clean after

\[ \small{ t = \frac{b - a}{v} \Leftrightarrow x = b } \]

• How long will it take a lake to decrease to \( 5\% \) of its current level of pollutant, assuming \( \, g(t) = 0 \) and \( P(x) \neq 0 \)?

• Lake has constant width and depth, so mass given by

\[ \small{ M=\iiint_E C(x, t) dV = A \int_a^b P(x-vt) dx } \]

• Find \( t \) such that

\[ \small{ \begin{aligned} \int_a^b P(x-vt) dx &= 0.05 \int_a^b P(x) dx \\ \int_{a-vt}^{b-vt} P(u)du &= 0.05 \int_a^b P(x) dx \\ \end{aligned} } \]

• How long will it take a lake to decrease to \( 5\% \) of its current level of pollutant, assuming \( \, g(t) = 0 \) and \( P(x) \neq 0 \)?
• Find \( t \) such that

\[ \small{ \int_{a-vt}^{b-vt} P(u)du = 0.05 \int_a^b P(x) dx } \]

• Left graph shows advection only case.
• Right graph shows diffusion only case.
  • Coloring shows level curves

Maple comands from book:

Maple comands from book:

• Math 366 Methods II: Ch21.6 (\( u_t = u_{xx} \), diffusion only).
• Grid with uniform spacing \( h \) along x-axis and \( \tau \) along t-axis.
• Forward time, centered space (FTCS) method

\[ \small{ \begin{aligned} \frac{\partial C}{\partial t} &= - v \frac{\partial C}{\partial x} \\ C_j^{n+1} & = C_j^n - \frac{v \tau }{2h}\left(C_{j+1}^n - C_{j-1}^n \right) \end{aligned}} \]

• Math 366 Methods II: Ch21.6 (Crank-Nicholson method).
• FTCS method is unstable; can use Lax method instead:

\[ \small{ \begin{aligned} C_j^{n+1} & = C_j^n - \frac{v \tau }{2h}\left(C_{j+1}^n - C_{j-1}^n \right) \,\,\, \mathrm{(FTCS)}\\ C_j^{n+1} & = \frac{1}{2}\left( C_{j+1}^n + C_{j-1}^n \right) - \frac{v \tau }{2h}\left(C_{j+1}^n - C_{j-1}^n \right) \,\,\, \mathrm{(Lax)} \end{aligned}} \]

• Method of Characteristics Video
  
  • B. Peters, UC Santa Barbara (2012)
• Method of Characteristics Honors Thesis
  • R. Malthaner, Baylor (2016)
• Finite Difference Methods for Advection Problems
  • P. Duffy, University College Dublin (2004)