In [3]:
import sys
sys.path.append('../')
import matplotlib.pyplot as plt
from porousmedialab.batch import Batch
import numpy as np
import seaborn as sns
%matplotlib inline
In [5]:
from porousmedialab.batch import Batch
bl = Batch(tend=40, dt=1)
bl.add_species(name='POC', init_conc=12e-3)
bl.add_species(name='CO2', init_conc=2e-3)
bl.add_species(name='Fe2', init_conc=0)
bl.add_species(name='CH4', init_conc=0)
bl.add_species(name='NO3', init_conc=1.5e-3)
bl.add_species(name='Fe3', init_conc=20e-3)
bl.add_species(name='SO4', init_conc=1.7e-3)
bl.constants['Km_NO3'] = 0.001e-3
bl.constants['Km_Fe3'] = 2e-3
bl.constants['Km_SO4'] = 0.3e-4
bl.constants['k1'] = 0.1
bl.rates['r_NO3'] = 'k1 * POC * NO3 / (Km_NO3 + NO3)'
bl.rates['r_Fe3'] = 'k1 * POC *Fe3 / (Km_Fe3 +Fe3) * Km_NO3 / (Km_NO3 + NO3)'
bl.rates['r_SO4'] = 'k1 * POC * SO4 / (Km_SO4 + SO4) * Km_Fe3 / (Km_Fe3 +Fe3) * Km_NO3 / (Km_NO3 + NO3)'
bl.rates['r_CH4'] = 'k1 * POC * Km_SO4 / (Km_SO4 + SO4) * Km_Fe3 / (Km_Fe3 +Fe3) * Km_NO3 / (Km_NO3 + NO3)'
bl.dcdt['POC'] = '- r_NO3 - r_Fe3 - r_SO4 - r_CH4'
bl.dcdt['NO3'] = '- 4 / 5 * r_NO3'
bl.dcdt['Fe3'] = '- 4 * r_Fe3'
bl.dcdt['Fe2'] = '4 * r_Fe3'
bl.dcdt['SO4'] = '- 1 / 2 * r_SO4'
bl.dcdt['CO2'] = 'r_NO3 + r_Fe3 + r_SO4 + 0.5 * r_CH4'
bl.dcdt['CH4'] = '1 / 2 * r_CH4'
bl.solve()
In [11]:
plt.rc('text', usetex=True)
SO4_t = [0.013, 1.013, 2.016, 3.005, 4.985, 6.851, 8.865, 11.882, 14.899, 18.713, 23.057, 28.923, 35.955]
SO4 = [1.643, 1.646, 1.565, 1.678, 1.484, 1.360, 1.121, 0.754, 0.336, 0.011, 0.014, 0.003, 0.002]
NO3_t = [0.036, 1.054, 1.988, 3.053, 5.008, 6.859, 9.018, 11.864, 14.857, 18.734, 23.065, 28.895, 35.898]
NO3 = [1.710, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Fe2_t = [0.017, 0.988, 2.003, 3.033, 4.987, 6.916, 8.909, 11.923, 14.896, 18.653, 23.070, 28.984, 35.975]
Fe2 = [1.588, 2.307, 5.167, 10.080, 15.267, 16.464, 18.682, 18.846, 16.893, 18.819, 19.503, 19.049, 19.254]
CO2_t = [0.156, 1.032, 2.070, 3.072, 5.069, 6.893, 8.919, 11.928, 14.948, 18.766, 23.078, 28.996, 35.988]
CO2 = [1.076, 2.109, 3.247, 5.008, 7.323, 8.342, 9.048, 9.648, 10.971, 10.958, 11.512, 10.448, 10.802]
CH4_t = [0.049, 1.016, 2.018, 3.011, 5.014, 6.872, 8.894, 11.903, 14.888, 18.724, 23.030, 28.954, 35.978]
CH4 = [0.005, 0.002, 0.005, 0.005, 0.005, 0.005, 0.022, 0.044, 0.077, 0.091, 0.168, 0.347, 0.637]
fig, ax1 = plt.subplots(figsize=(5,3), dpi=300)
ax2 = ax1.twinx()
ax1.plot(bl.time, 1e+3*bl.NO3.concentration[0], c=sns.color_palette("deep", 10)[0], lw=3, label='$NO_3$')
ax1.scatter(NO3_t, NO3, c=sns.color_palette("deep", 10)[0], lw=1)
ax1.plot(bl.time, 1e+3*bl.SO4.concentration[0], c=sns.color_palette("deep", 10)[1], lw=3, label='$SO_4$')
ax1.scatter(SO4_t, SO4, c=sns.color_palette("deep", 10)[1], lw=1,)
ax1.plot(bl.time, 1e+3*bl.CH4.concentration[0], c=sns.color_palette("deep", 10)[2], lw=3, label='$CH_4$')
ax1.scatter(CH4_t, CH4, c=sns.color_palette("deep", 10)[2], lw=1,)
ax2.plot(bl.time, 1e+3*bl.Fe2.concentration[0], c=sns.color_palette("deep", 10)[3], lw=3, label='$Fe(II)$')
ax2.scatter(Fe2_t, Fe2, c=sns.color_palette("deep", 10)[3], lw=1,)
ax2.plot(bl.time, 1e+3*bl.CO2.concentration[0], c=sns.color_palette("deep", 10)[4], lw=3, label='$\sum CO_2$')
ax2.scatter(CO2_t, CO2, c=sns.color_palette("deep", 10)[4], lw=1,)
ax2.grid(False)
ax1.grid(lw=0.2)
ax1.set_ylim(0, 2.5)
ax1.set_xlim(0, 40)
ax2.set_ylim(0, 25)
ax1.set_ylabel('$NO_3$, $SO_4$, $CH_4$, [mM]')
ax2.set_ylabel('$Fe(II)$, $\sum CO_2$, [mM]')
ax1.set_xlabel('Time, [days]')
ax1.legend(frameon=1, loc=2)
ax2.legend(frameon=1, loc=1)
Out[11]:
In [ ]: