AN_tmf_assignment_pha.ipynb

Kernel: Python 3 (Ubuntu Linux)

Table of Contents

• 1  Assignment
  • 1.1  Objectives and description**
  • 1.2  Problem description and conceptual model
  • 1.3  Identification of sulfate sources
    • 1.3.1  Parameters
  • 1.4  Evolution of the mass of water
  • 1.5  Evolution of the mass of sulfates
  • 1.6  Mass balance
    • 1.6.1  Associated ODE
    • 1.6.2  Influence of the timestep

Assignment

Objectives and description**

The objective of this assignment is to create a python program that computes the concentration of sulfate in time in the water in a tailings management facility (TMF) at a mine site.

The goal is to describe the computational approach to the resolution of a typical "0D" problem. In the end, the goal is to understand how the computational approach corresponds to the resolution of an ODE. We will investigate that ODE and study its solutions. That will allow us to compare the accuracy of the numerical solution, because mathematics allowed us to know, in this particular case, the exact solution of a problem.

This assignment is built so that the concepts of the accuracy of a solution are naturally introduced. In the end, you will realize that you will have implemented what is called the forward Euler approach.

Problem description and conceptual model

Sulfate $SO_{4}^{2 -}$ is one of the “major” dissolved constituents in most terrestrial waters (the others being sodium, calcium, magnesium, bicarbonate and chloride). Sulfate is not poisonous, but much in the same way that too much salt makes water undrinkable, too much sulfate can be too and is thus an environmental concern.

Furthermore, high sulfate concentrations are usually associated with acidic conditions (think about sulfuric acid producing sulfates). In acidic conditions, heavy metals (which are poisonous) might be released from rock dissolution, ... Therefore, monitoring and predicting the sulfate content in groundwater is an essential environmental concern.

We will go back to the TMF problem, where we want to know and assess the sulfate concentration evolution through time.

Identification of sulfate sources

To compute the concentration of sulfate in the TMF water, we need to know where it is coming from and going to and use that knowledge to develop a conceptual model describing the fate of sulfate in our system. That conceptual model will be translated into a computational model.

The figure below is a conceptual sketch of the main sources and losses of sulfate in the TMF water:

1. From the pit. To keep the pit dry, water is constantly pumped out of the pit and put into the TMF at a rate $Q_{\text{pit}}$ (L/s of water). The pit water has a sulfate concentration of $c_{\text{pit}}$ in mg/L.

2. From the mill. The tailings are pumped into the TMF as a slurry of tailings particles (fine, sand-size ground rock) and water. Assume that the water from the mill enters the TMF at a rate of $Q_{\text{mill}}\$ (L/s of water), with a sulfate concentration of $c_{\text{mill}}$ (mg/L).

3. Diffusing from the tailings porewater at the bottom of the TMF into the water column above. Sulfate dissolves from the tailings particles into the adjacent porewater by oxidation of the sulfide minerals in the particles. Because the ratio of rock to water is high in the tailings sediments at the bottom, the porewater sulfate concentration in the tailings, $c_{\text{pore}}$, is always higher than in the water in the TMF. Accordingly, sulfate tends to diffuse into the TMF water from the bottom porewater at a rate proportional to the difference in concentration between the porewater concentration and the concentration in the TMF, $c_{\text{TMF}}.$ That is, the flux of sulfate from the bottom can be written (a positive quantity when sulfate is entering into the TMF):

$j_{\text{pore}} = k\left( c_{\text{pore}} - c_{\text{TMF}} \right)$

where $j_{\text{pore}}\ mg/(s \cdot m^{2})$ is the flux rate of sulfate per unit area of the bottom of the TMF and $k$ is the flux coefficient with units of $L/(s \cdot m^{2})$ . To compute $J_{\text{pore}}$, the total mass flux rate for the whole TMF, with units of $\frac{\text{mg}}{s}$, we must multiply the rate per unit area $j_{\text{area}}$ by the total area of the TMF bottom $A_{\text{bottom}}\$ in $m^{2}$:

$J_{\text{pore}} = A_{\text{bottom}}{\times j}_{\text{pore}}$

4. Leaving via the discharge ditch. Sulfate leaves the TMF with the water that is discharged at a rate $Q_{\text{discharge}}$ (L/s) from the TMF to the environment. SIMPLIFYING ASSUMPTION: if we assume that the water in the TMF is well mixed at all times, then the concentration in the TMF, $c_{\text{TMF}}$ is the same at all points in the TMF, and therefore the water leaving the TMF has that same concentration $c_{\text{TMF}}\$.

Parameters

You will need these parameters for your model

Symbol	Symbol Units	Description	Value
$c_{\text{pit}}$	$mg/L$	Concentration of sulfate in pit water (assume constant)	$50$
$Q_{\text{pit}}$	$L/s$	Flow rate of water from the pit into TMF (assume constant)	30
$c_{\text{mill}}$	$mg/L$	Concentration of sulfate in mill water (assume constant)	$700$
$Q_{\text{mill}}$	$L/s$	Flow rate of water from mill into TMF (assume constant)	14
$Q_{\text{discharge}}$	$L/s$	Flow rate of water from TMF to environment (assume constant)	44
$c_{\text{pore}}$	$mg/L$	Concentration of sulfate in porewater at bottom of pond	2000
$k$	$L/(s \cdot m^{2})$	Flux coefficient from porewater to water column	$2.5 \times 10^{- 5}$
$A_{\text{bottom}}$	$m^{2}$	Total area of TMF bottom	$3 \times 10^{5}$
$V_{\text{TMF}}$	$L$	Volume of water in TMF at start of simulation	$8.1 \times 10^{6}$
$c_{0}$	$mg/L$	Concentration of sulfate in TMF water at start of simulation	$93$

Evolution of the mass of water

First, we will check the evolution of the volume of water in the TMF. Let us denote by $V(t)$ the volume of water in the TMF at every time step.

The flux from the pit and the mill are positive source terms for the volume of water, while the discharge is a negative source term (a sink term).

Over a certain time $\Delta t$, the amount of water (in L) which is coming in/out the TMF are denoted $S_{\text{pit}}$, $S_{\text{mill}}$, $S_{\text{dis}}$ and are equal to

$$\begin{equation} \left\lbrace \begin{array}{lll} S_{\text{pit}} & = & Q_{\text{pit}} \Delta t \\ S_{\text{mill}} & = & Q_{\text{mill}} \Delta t \\ S_{\text{dis}} & = & Q_{\text{dis}} \Delta t \\ \end{array} \right. \end{equation}$$

Therefore, the volume $V$ of water in the TMF over the period [$t$;$t+\Delta t]$ is:

$$\begin{equation} \begin{array}{llll} & V(t+\Delta t) & = & V(t) + S_{\text{pit}} + S_{\text{mill}} - S_{\text{dis}} \\ \Longleftrightarrow & V(t+\Delta t) & = & Q_{\text{pit}} \Delta t + Q_{\text{mill}} \Delta t - Q_{\text{dis}} \Delta t \\ \Longleftrightarrow & \frac{V(t+\Delta t)-V(t)}{\Delta t} & = & Q_{\text{pit}} + Q_{\text{mill}} - Q_{\text{dis}} \\ \end{array} \end{equation}$$

How does the volume of water change with time?

from numpy.testing import assert_almost_equal
Q_pit = 30
Q_mill = 14
Q_dis = 44

V_change = Q_pit+Q_mill-Q_dis
print(V_change)

0

Considering the values given, we have:

$$\begin{equation} \frac{V(t+\Delta t)-V(t)}{\Delta t} = 30 + 14 - 44 = 0 \Longleftrightarrow V(t) = V_0 \end{equation}$$

Indicating that the volume of water is constant through time!

Let's assigne a variable whose value is the volume of water in the TMF:

V0 = 0
#  Change the value of V0 to match the given data

V0=8.1E9

assert_almost_equal(V0,8.1e9,decimal=1)

Evolution of the mass of sulfates

We have 4 main processes impacting the quantity of sulfates in the TMF. We will compute the evolution of this mass over time. First, initialize the required variable (initial mass, and the different input parameters).

c0 = 93  # initial concentration

c_pit = 50
c_pore = 2000
c_mill = 700

k = 2.5e-5
Area = 3e5

What is the initial mass of sulfates? (in mg)

m0 = 0
# Change the value of m0 so that it corresponds to the initial mass (in mg) of sulfates in the TMF.

m0=c0*V0
print(m0)

753300000000.0

assert_almost_equal(m0,7.533e11,decimal=1)

1. Advective flux from the pit and the mill

We have a source of sulfates from the pit. Over a certain time $\Delta t$, the volume of water coming from the pit is $S_{\text{pit}} = Q_{\text{pit}} \Delta t$. This volume of water has a sulfate concentration corresponding to the pit concentration. So, the total mass of sulfate from the pit which arrived in the pit over the same period, $M_{\text{pit}}$ is simply: $$\begin{equation} M_{\text{pit}} = S_{\text{pit}} c_{\text{pit}} = Q_{\text{pit}} c_{\text{pit}} \Delta t \end{equation}$$

You can always check the units: $$\begin{equation} \overbrace{M_{\text{pit}}}^{mg} = \overbrace{Q_{\text{pit}}}^{L/s} \, \overbrace{c_{\text{pit}}}^{mg/L} \, \overbrace{\Delta t}^{s} \end{equation}$$

Compute the mass of sulfates after 1 day if the pit is the only source of sulfate. Create a variable S_pit corresponding to the mass (mg) of sulfates brought by the pit in one day.

dt = 1  # day
Seconds_in_a_day = 24 * 3600
S_pit = 0

# YOUR CODE HERE
S_pit = dt*Seconds_in_a_day*c_pit*Q_pit
print(S_pit)

129600000

assert_almost_equal(S_pit,129600000,decimal=1)

2. Advective flux from the mill

If we only consider the mill as a sulfate source, how will the mass of sulfates evolve in one day. Assign this value in variable named S_mill (in mg).

S_mill = 0

S_mill=dt*Seconds_in_a_day*Q_mill*c_mill
print(S_mill)

846720000

#assert that S_mill = 846720000
assert_almost_equal(S_mill,846720000,decimal=1)

3. Discharge flux from the TMF

The same development can be performed for the discharge of water. The volume of water leaving the TMF over a certain period $\Delta t$ is $S_{\text{dis}} = Q_{\text{dis}} \Delta t$. The concentration of sulfates in this volume corresponds to the concentration of sulfates in the TMF $c_{\text{TMF}}$, so that we can write:

$$\begin{equation} M_{\text{dis}} = Q_{\text{dis}} c_{\text{TMF}} \Delta t \end{equation}$$

How will the discharge change the mass of sulfates over one day? Create a variable called S_dis which stores that value (in mg).

S_dis = 0

# YOUR CODE HERE
S_dis = dt*Seconds_in_a_day*c0*Q_dis
print(S_dis)

353548800

assert_almost_equal(S_dis,353548800,decimal=1)

4. First order mass transfer

The mass flux was written as: $$\begin{equation} J_{\text{pore}} = \overbrace{A_{\text{bottom}}}^{m^2} \overbrace{k}^{L . s^{-1} . m^{-2}} \left( \overbrace{c_{\text{pore}} - c_{\text{TMF}}}^{mg . L^{-1}} \right) \end{equation}$$ which is in the units of mg/s. If this flux is constant over time, the added mass of sulfates over a certain time $\Delta t$ corresponds to the latter quantity multiplied by the the quantity $J_{\text{pore}}$ multiplied by $\Delta t$:

$$\begin{equation} M_{\text{pore}} = A_{\text{bottom}} k \left( c_{\text{pore}} - c_{\text{TMF}} \right) \Delta t \end{equation}$$

How will this mass transfer impact the mass of sulfate? Create a variable S_pore which stores the value (in mg) of the sulfate brought by this first order mass transfer.

S_pore = 0

# YOUR CODE HERE
S_pore=Area*k*(c_pore-c0)*dt*Seconds_in_a_day
print(S_pore)

1235736000.0

#ASSERT THAT S_PORE = 1235736000.0
assert_almost_equal(S_pore,1235736000,decimal=1)

Now, compute the total evolution of the mass of sulfates after one day, and print its value and how it has changed. Store that value (in mg) in a variable called m1.

m1 = 0

m1=(m0+S_mill+S_pit-S_dis+S_pore)*dt
print(m1)

755158507200.0

change_in_mass = m1-m0
print(change_in_mass)

1858507200.0

assert_almost_equal(m1,755158507200,decimal=1)

This is the mass after 1 day (in mg). What is the concentration after one day? Store that value in a variable called c1 (in mg/L).

c1 = 0

c1=m1/V0
print(c1)

93.22944533333333

#assert that c1=93.2294453333
assert_almost_equal(c1,93.2294453333,decimal=3)

We can see that concentration has increased by only 0.2%.

We have computed the evolution of the mass of sulfates in the TMF after one day. How would you do that after another day?

We will use a loop to do that calculation, to compute the evolution of the concentration over a large period of time (10 years, for example), using daily timesteps again.

At the end, we want to plot the evolution of the concentration, so we have to store its values in an array, whose size corresponds to the amount of timesteps required.

Using daily timesteps, what size would that array need to be if we want to monitor the concentration for 10 years. That size is an integer called n.

Tf = 10 #years

Days_in_a_year = 365
n=Tf*Days_in_a_year+1
print(n)

3651

assert_almost_equal(n,3651,decimal=1)

import matplotlib.pyplot as plt
import numpy as np

# Let us assign some arrays to describe our problem
m = np.zeros(n, float) # represents the mass of sulfates (mg) at each time
c = np.zeros(n, float) # represents the concentration of sultates (mg/L) at each time
Spit = np.zeros(n, float) # represents at each time, the amount of sulfates (in kg) brought by the pit
Smill = np.zeros(n, float) # represents at each time, the amount of sulfates (in kg) brought by the mill
Sdis = np.zeros(n, float) # represents at each time, the amount of sulfates (in kg) discharged
Spore = np.zeros(n, float) # represents at each time, the amount of sulfates (in kg) brought by the tailing
time = np.zeros(n, float) # represents the time at each timestep
discharge = np.zeros(n, float) # represents the cumulative discharge of sulfates (in kg)

# First, we initialize the initial mass and concentrations.

t0=0
mill0=S_mill/V0
pit0=S_pit/V0

m[0] = m0
c[0] = c0
time[0]=t0
Smill[0]=mill0
Spit[0]=pit0

print(pit0)
type(pit0)

0.016

float

discharge[0]=S_dis/1000000

Then, we need to compute these terms at every timesteps, by looping over the different times

for i in range(1, n):
    Spore[i]=Area*k*(c_pore-c[i-1])*86400/1000000
    m[i]=m[i-1]+Spore[i]*1000000+S_mill+S_pit-(Sdis[i-1]*1000000)
    c[i]=m[i]/V0
    Sdis[i]=(c[i-1]*Q_dis*86400)/1000000
    time[i]=time[i-1]+1
    Smill[i]=Smill[i-1]
    Spit[i]=Spit[i-1]
    discharge[i]=discharge[i-1]+Sdis[i]
      

print(c[i-1])

454.4974144339691

print(discharge[n-1])

4586022.146993127

assert_almost_equal(c[n-1],454.46947737053875,decimal=0)

assert_almost_equal(discharge[n-1],4.584599422624083e+06,decimal=1)

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-42-8161f5926ad9> in <module>()
----> 1 assert_almost_equal(discharge[n-1],4.584599422624083e+06,decimal=1)

/usr/local/lib/python3.6/dist-packages/numpy/testing/_private/utils.py in assert_almost_equal(actual, desired, decimal, err_msg, verbose)
    582         pass
    583     if abs(desired - actual) >= 1.5 * 10.0**(-decimal):
--> 584         raise AssertionError(_build_err_msg())
    585 
    586 
AssertionError: 
Arrays are not almost equal to 1 decimals
 ACTUAL: 4586022.146993127
 DESIRED: 4584599.422624083

import numpy as np
from matplotlib import pyplot as plt

fig, (ax1) = plt.subplots(1, 1, figsize=(8, 4))
ax1.plot(time, c, label="Pond Concentration")
ax1.set(xlabel="Time (days)", ylabel="Concentration (mg/L)")
ax1.legend()

fig, (ax1) = plt.subplots(1, 1, figsize=(8, 4))
ax1.plot(time, Smill, label="mill contribution")
ax1.plot(time, Spit, label="pit contribution")
ax1.plot(time, Spore*1000000/V0, label="pore contribution")
ax1.plot(time, -Sdis*1000000/V0, label="discharge contribution")
#ax1.plot(time, c, label="Pond Concentration")
ax1.set(xlabel="Time (days)", ylabel="Concentration (mg/L)")
ax1.legend()

The calculation is now performed, let us look at the results using matplotlib. Design 3 vertically stacked plots

1. The first should show the evolution of the concentration over time

2. The second should show the relative importance of the different sources/sinks

3. The third should show the cumulative discharge (in kg) of sulfates

ig, (ax1) = plt.subplots(1, 1, figsize=(8, 4))
ax1.plot(time, discharge, label="Cumulative Discharge")
ax1.set(xlabel="Time (days)", ylabel="Concentration (mg/L)")
ax1.legend()

Mass balance

Associated ODE

We can (and you will!) show that this problem can be described by a linear non-homogeneous 1st order ODE. We can (and you will) verify that the solution to this ODE is given by:

$$\begin{equation} c_{\text{TMF}}(t) = \frac{Q}{\lambda} + A \mathrm{exp}(-\lambda t) \end{equation}$$

where $$\begin{equation} \left\lbrace \begin{array}{lll} Q & = & \frac{Q_{\text{pit}}c_{\text{pit}} + Q_{\text{mill}}c_{\text{mill}}+ A_{\text{bottom}}kc_{\text{p}}}{V_0} \\ \lambda & = & \frac{A_{\text{bottom}}k+Q_{\text{dis}}}{V_0} \end{array} \right. \end{equation}$$

Find the value of A so that the initial condition is satisfied.

% A = 0 Assign A to its real value so that the previous solution matches the initial ccondition of the problem

print(k)
Area=3E5
print(Area)
print(Q_dis)
print(Area)
print(V0)
print(Q_dis)

lam = (Area * k + Q_dis) / V0
print(lam)

Q = (Q_pit*c_pit+Q_mill*c_mill+Area*k*c_pore)/V0
print(Q_pit)
print(c_pit)
print(Q_mill)
print(c_mill)
print(c_pore)
print(Q)

# Parameters definition
import numpy as np
from matplotlib import pyplot as plt

c_pit = 50
c_pore = 2000
c_mill = 700
Q_pit = 30
Q_mill = 14
Q_dis = 44
V0 = 8.1e9
c0 = 93  # initial concentration
k = 2.5e-5
Area = 3e5
#
# Eqn 13
#
c_inf = (Q_pit * c_pit + Q_mill * c_mill + Area * k * c_pore) / (Q_dis + Area * k)
A = c0 - c_inf
lam = (Area * k + Q_dis) / V0  # units are here in s^{-1}
#
# Eqn 15
#
Q = (Q_pit * c_pit + Q_mill * c_mill + Area * k * c_pore) / V0
print(lam)
print(Q)
print(A)

c0=93
print(c0)
A = c0-(Q/lam)
print(A)

Again, we want a triple vertical plots (with time as x-axis).

1. The first plot shows the computed solution vs the analytical solution

2. The second plot shows the absolute error (computed_solution - real solution)

3. The third plot shows the relative error (computed_solution-real_solution)/real_solution

To do that, you need to compute the real solution and store it in an array. You also have to compute the error and relative error. The next cell should contain the analytical calculation of the real solution and every plot material.

# This is the routine to plot the exact solution for the concentration at every day
dt = 1  # day
Tf = 10  # years
n = int(1 + 365 * Tf / dt)  # int() is to convert thhe result into an integer
c_real = np.zeros(
    n, float
)  # represents the concentration of sultates (mg/L) at each time
timeY = np.zeros(n, float)  # time in years
c_real[0] = c0
c_asymptotic = c_inf * np.ones(
    n, float
)  # this is an array of size n full, whose values are all c_inf
Forward_Euler = np.zeros(n, float)
Backward_Euler = np.zeros(n, float)
Abs_err = np.zeros(n, float)
Rel_err = np.zeros(n, float)
Forward_Euler[0] = c0
Backward_Euler[0] = c0


for i in range(n - 1):
    time[i + 1] = (i + 1) * dt / 365
    c_real[i + 1] = c_inf + (c0 - c_inf) * np.exp(-lam * (i + 1) * dt * 3600 * 24)
    Forward_Euler[i + 1] = Forward_Euler[i]+((Q-(lam*Forward_Euler[i]))*(dt*86400))
    Backward_Euler[i + 1] = (Backward_Euler[i]+(Q*dt*24*3600))/(1+lam*dt*86400)
    Abs_err[i + 1] = c[i+1] - c_real[i + 1]
    Rel_err[i + 1] = (c[i + 1] - c_real[i + 1]) / c_real[i + 1]

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 8))
ax1.plot(time, c, label="Computed Concentration")
#ax1.plot(time, Forward_Euler, label="Concentration - Analytical")
ax1.plot(time, c_real, label="Concentration - C_Real")
ax1.set(xlabel="Time (days)", ylabel="Concentration (mg/L)")
ax1.legend()

ax2.plot(time, 100 * Abs_err, label="Error - Absolute")
#ax2.plot(time, 100 * Rel_err, label="Error - Relative")
ax2.set(xlabel="Time (years)", ylabel="Error (%)")
ax2.legend()

#ax3.plot(time, 100 * Abs_err, label="Error - Absolute")
ax3.plot(time, 100 * Rel_err, label="Error - Relative")
ax3.set(xlabel="Time (years)", ylabel="Error (%)")
ax3.legend()

Influence of the timestep

The error is arising from the timestepping approach we have used. During one timestep, we have computed the discharge flux and the source-term from the "pore" as if the concentration $c_{\text{TMF}}$ was constant through the timestep. This is actually not true. So, the bigger the timesteps, the bigger the error.

Try the same methods used above for different timesteps (0.1 day, 1 day, 10 days, 50 days) and compare each of these solutions to the real solution and comment on the error.

# YOUR CODE HERE
raise NotImplementedError()

How would you change the method if the volume of water was not constant?

# YOUR CODE HERE
raise NotImplementedError()