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.
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.
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:
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.
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).
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}}$
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$ |
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
Therefore, the volume $V$ of water in the TMF over the period [$t$;$t+\Delta t]$ is:
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)
Considering the values given, we have:
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)
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)
assert_almost_equal(m0,7.533e11,decimal=1)
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:
You can always check the units:
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)
assert_almost_equal(S_pit,129600000,decimal=1)
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)
#assert that S_mill = 846720000 assert_almost_equal(S_mill,846720000,decimal=1)
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:
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)
assert_almost_equal(S_dis,353548800,decimal=1)
The mass flux was written as: 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$:
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)
#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)
change_in_mass = m1-m0 print(change_in_mass)
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)
#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)
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)
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])
print(discharge[n-1])
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
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()
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:
where
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).
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()
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()