o POSSIBLE SOLUTION(S)

+ Problem 1. In and Out the Same --- Instantaneous Mixing

- Key Assumptions:
Pollutants dilute and mix uniformly throughout lake instantaneously!
No evaporation (and thus, the lake volume is constant, V)

Then, if
P(t) = pounds of pollutants at time t (in minutes),
V = volume of the lake (in gallons),
Po = initial amount of pounds of pollutants at time t (in minutes);

then the differential equation describing the amount of polution is:


dP/dt = (incoming rate) - (outgoing rate)
= 6000 - (P(t) [lb])/(V [gal])(5000 [gal/min] )
= 6000 - (5000/V) P(t)

P(0) = Po.


Note that since the right hand side of the DE is time dependent (i.e., the DE is autonomous),

P_eq = 6000 1000000/5000 = 1.2 million pounds

is an equilibrium solution. Furthermore, it's attracting, i.e., as t -> infinity,
P(t) -> P_eq.

- Here's how to solve the IVP (assuming the lake starts out with fresh water...) 

Input := 

Clear[p, t];

Input := 

V = 1000000;
eqn = p'[t] == 6000 - 5000/V p[t];
ic = p[0] == 0;
ans = Flatten[DSolve[{eqn, ic}, p[t], t]];

Input := 

sol = p[t] /. ans;

Input := 

psol = Plot[sol, {t, 0, 100}]

Output =

-Graphics-

+ Problem 2. Two layers.

- Assume that
p1[t] = pounds of pollutant in the top layer;
p2[t] = pounds of pollutant in the bottom layer;
d12 = rate (in gallons/minute) that pollutants transfer from top to bottom;
d21 = rate (in gallons/minute) that pollutants transfer from bottom to top;
V = volume of lake (in gallons);

- Then, the differential equations that govern the pollutants in each layer are: 

Input := 

V = 1000000;
e1 = p1'[t] == 6000 - 5000/(V/2) p1[t] - d12/(V/2) p1[t] +
                 d21/(V/2) p2[t];
e2 = p2'[t] == d12/(V/2) p1[t] - d21/(V/2) p2[t];

- The Fixed Point Analysis

- shows that the unique fixed solution 

Input := 

eq1 = 6000 - 5000/(V/2) p1[t] - d12/(V/2) p1[t] +
            d21/(V/2) p2[t]==0;
eq2 = d12/(V/2) p1[t] - d21/(V/2) p2[t]==0;
Flatten[Solve[{eq1, eq2}, {p1[t], p2[t]}]]

Output =

          600000 d12
{p2[t] -> ----------, p1[t] -> 600000}
             d21

- is attracting. Note that the limiting weight of pollutants in the top layer (p1[t]) is independent of the diffusivities for this fixed solution.

- For those familiar with matrices, one can show that the solutions are attracting by use of the Jacobian matrix which is the matrix associated with the linear system: 

Input := 

jac = {{- 5000/(V/2) - d12/(V/2), d21/(V/2)},
 {d12/(V/2), -d21/(V/2)}}

Output =

     1      d12     d21       d12     -d21
{{-(---) - ------, ------}, {------, ------}}
    100    500000  500000    500000  500000

- and the eigenvalues of this Jacobian matrix are 

Input := 

Simplify[Eigenvalues[jac]]

Output =

{(-500000 - 100 d12 - 100 d21 - 
 
                                                   2
     Sqrt[-200000000 d21 + 10000 (5000 + d12 + d21) ]) / 100000000, 
 
  (-500000 - 100 d12 - 100 d21 + 
 
                                                   2
     Sqrt[-200000000 d21 + 10000 (5000 + d12 + d21) ]) / 100000000}

- which always have negative real part. So, the equilibrium solutions are stable (attracting.)

- Initial Conditions and Corresponding Solutions

If we assume that the lake is initially clean, we have the following initial conditions: 

Input := 

ic1 = p1[0] == 0;
ic2 = p2[0] == 0;

Here's the general solutions:
sol1 and sol2 (in terms of the unspecified parameters d12, d21, and V.) 

Input := 

ans = Flatten[DSolve[{e1, e2, ic1, ic2}, {p1[t], p2[t]}, t]];

Input := 

sol1 = p1[t] /. ans[[1]];
sol2 = p2[t] /. ans[[2]];

- Here are some sample graphs of solutions and comments on the results...

Here we show that when d12 = d21 is small, p1(t) and p2(t) are quite different and p1(t) > p2(t) (which makes sense since most particulate matter remains near the surface where it came into the lake water), while as d12 = d21 -> infinity, mixing becomes instantaneous throughout the lake and p1(t) and p2(t) coalesce. The animation further demonstrates that p2(t) increases from nearly zero (the initial weight of pollutants) as the diffusivity increases, while p1(t) decreases with increasing diffusivity since more of the pollutants are able to cross to the other layer.

(In the animation, the parameter d12 = d21 increases exponentially. The thicker curve is p1(t).) 

Input := 

Do[Plot[{sol1 /. {d12 -> N[Exp[i]], d21 -> N[Exp[i]]},
      sol2 /. {d12 -> N[Exp[i]], d21 -> N[Exp[i]]}},
 {t, 0, 100},
 PlotStyle -> {AbsoluteThickness[2],
               AbsoluteThickness[1]},
 AxesLabel -> {"minutes", "pollutant weight"},
 PlotRange -> {{0, 100}, {0, 350000}}],
    {i, 5, 13}]

- Here, we consider the affects of changing the ratio of the diffusivities. This is done by letting alpha = d12/d21 and then fixing -- say d21 = 1 while varying d12 both above and below 1000.

In the animation, the parameter alpha increases exponentially. The thicker curve is p1(t).

Note that when alpha is small, most pollutants remain in the top layer, but when alpha is very large, most pollutants end up in the bottom layer. 

Input := 

Do[Plot[{sol1 /. {d21 -> 1., d12 -> N[Exp[alpha]]},
      sol2 /. {d21 -> 1., d12 -> N[Exp[alpha]]}},
 {t, 0, 100},
 PlotStyle -> {AbsoluteThickness[2],
               AbsoluteThickness[1]},
 AxesLabel -> {"minutes", "pollutant weight"},
 PlotRange -> {{0, 100}, {0, 350000}}],
    {alpha, 6, 12, 0.5}]

- Below we compare the graphs of the solutions when the lake is assumed to be in two layers vs. when the lake is assumed to mix entirely instantaneously. The bold graph corresponds to no layers (instantaneous mixing throughout), while the thinner graph corresponds to instantaneous mixing within each layer.

This animation shows that as the diffusivities increase and thus mixing occurs more readily, that the simplified model becomes a reasonable approximation. The parameters d12 = d21 increase exponentially. This is an important result for the company dumping pollutants since the simplified model over estimates the amount of damage done to the environment.

The solution of the simplified model (from #1) is the thicker curve and does not change. 

Input := 

Do[Plot[{sol1+sol2 /. {d21 -> N[Exp[i]], d12 -> N[Exp[i]]},
      sol},
 {t, 0, 100},
 PlotStyle -> {AbsoluteThickness[1],
               AbsoluteThickness[2]},
 AxesLabel -> {"minutes", "pollutant weight"},
 PlotRange -> {{0, 100}, {0, 500000}}],
    {i, 5, 12, 0.5}]

+ Problem 3. Multiple Layers and Cells.

- There are many levels of sophistication that one can impose here. To start simply, one might partition the lake by making lots of horizontal slabs (say, of equal volume) and assuming that there is instantaneous mixing occurs within each horizontal slab. If

p_i(t) = pounds of pollutants in slab i,

then the model would be

p_1'(t) = 6000 - (n/V) { 5000 p_1(t) + d21 p_2(t) - d12 p_1(t) }
p_2'(t) = (n/V) { d12 p_1(t) + d32 p_3(t) - (d21 + d23) p_2(t) }
p_3'(t) = (n/V) { d23 p_2(t) + d43 p_4(t) - (d32 + d34) p_3(t) }
p_4'(t) = (n/V) { d34 p_3(t) + d54 p_5(t) - (d43 + d45) p_4(t) }
.
.
.
p_n'(t) = ...

It's a bit easier to see from the Jacobian matrix that this system is tridiagonal.

At a higher level of complexity, one can extend to two or even three dimensions. In these cases, the modeler must decide how pollutants can diffuse in higher dimensions. A reasonable assumption is that there is only diffusion across borders of dimension one less than the dimension of the cell.

For example, in two dimensions (the lake is partitioned into thin cylinders with small rectangular cross sections), pollutants move across the lines (1D) between neighboring cells, but not vertices at the diagonals.

- Code for Grid 

Input := 

Show[Graphics[{Line[{{0, 0}, {4, 0}}],
               Line[{{0, 0}, {0, 4}}],
               Line[{{0, 4}, {4, 4}}],
               Line[{{4, 4}, {4, 0}}],
               Line[{{1, 0}, {1, 4}}],
               Line[{{2, 0}, {2, 4}}],
               Line[{{3, 0}, {3, 4}}],
               Line[{{0, 1}, {4, 1}}],
               Line[{{0, 2}, {4, 2}}],
               Line[{{0, 3}, {4, 3}}],
               Text["Diffusion This Way", {1.3, 1.4}],
               Line[{{1.8, 1.5}, {2.2, 1.5}}],
               Line[{{2.1, 1.4}, {2.2, 1.5}}],
               Line[{{2.1, 1.6}, {2.2, 1.5}}],
               Text["But Not This Way", {2.3, 2.4}],
               Line[{{2.2, 2.2}, {1.8, 1.8}}],
               Line[{{1.8, 1.8}, {1.85, 1.9}}],
               Line[{{1.8, 1.8}, {1.9, 1.85}}]},
              AspectRatio -> Automatic]]