In [ ]:
%matplotlib inline
First create a new file by "Saving-as" NPDZ.py witha new name... lets call the new file NPZD2_Mussels.py
From now on, we are going to copy-paste from Mussels.py to the new NPZD2_Mussels.py
Copy-paste parameters from Mussels.py to NPZD2_Mussels.py Copy-paste the "Mussel" parameters just below the "NPZD2" parameters. The new code should look like this:
In [ ]:
# Parameters
par = {}
par['mu0'] = 0.69
par['kNO3'] = 0.5
par['kNH4'] = 0.5
par['alpha'] = 0.125
par['gmax'] = 0.6 #Original 0.6
par['kP'] = 0.44
par['mP'] = 0.15
par['tau'] = 0.005
par['thetaMax'] = 0.053
par['beta'] = 0.75
par['lBM'] = 0.1
par['lE'] = 0.1
par['mZ'] = 0.25
par['rSD'] = 0.3 # Original 0.03
par['rLD'] = 0.1 # # Original 0.01
par['nmax'] = 0.05
par['kI'] = 0.1
par['I0'] = 0.0095
# Parameters Mussels
par['AE_P'] = 0.9
par['AE_D'] = 0.2
par['AE_Z'] = 0.3
par['Bpub'] = 0.43
par['Fmax_ref'] = 0.025
par['GT'] = 0.44
par['KTempH'] = 0.1
par['KTempL'] = 0.5
par['KSaltL'] = 0.25
par['KOxyL'] = 0.02
par['KFood'] = 1.
par['KRE'] = 0.86
par['OxyL'] = 17.5
par['Rm'] = 0.002
par['SaltL'] = 10.
par['TempH'] = 25.
par['TempL'] = -4.
par['beta'] = 0.12
par['epsilonP'] = 1.
par['epsilonD'] = 0.5
par['epsilonZ'] = 0.3
Copy-paste InitCond from Mussels.py to NPZD2_Mussels.py Copy-paste the "Mussel" InitCond just below the "NPZD2" InitCond.
Note that now there are duplicate InitCond... do not copy the Mussel InitCond that are already in NPZD2 (i.e. Temp, Phy, Zoo, SDet)
The new code should look like this:
In [ ]:
# Initial conditions
InitCond = {}
InitCond['Phy'] = 0.2
InitCond['Zoo'] = 0.1
InitCond['SDet'] = 1.
InitCond['LDet'] = 1.
InitCond['NH4'] = 0.1
InitCond['NO3'] = 7.
InitCond['Temp'] = 6.
#InitCond['O2'] = 0.5
InitCond['Soma'] = 0.01
InitCond['Gonad'] = 0.
InitCond['Salt'] = 30. #Salinity
InitCond['Oxy'] = 30. #Oxygen
Copy-paste zero vectors from Mussels.py to NPZD2_Mussels.py Copy-paste the "Mussel" zero-vectors just below the "NPZD2" zero-vectors.
Same as above, note that now there are duplicate zero-vectors... do not copy the Mussel zero-vectors that are already in NPZD2 (i.e. Temp, Phy, Zoo, SDet)
The new code should look like this:
In [ ]:
# Create zero-vectors
Phy = np.zeros((NoSTEPS,),float) # makes a vector array of zeros (size: NoSTEPS rows by ONE column)
Zoo = np.zeros((NoSTEPS,),float) # same as above
SDet = np.zeros((NoSTEPS,),float) # Biomass - same as above
LDet = np.zeros((NoSTEPS,),float) # same as above
NH4 = np.zeros((NoSTEPS,),float) # same as above
NO3 = np.zeros((NoSTEPS,),float) # same as above
I = np.zeros((NoSTEPS,),float) # same as above
mu = np.zeros((NoSTEPS,),float) # same as above
f_I = np.zeros((NoSTEPS,),float) # same as above
L_NO3 = np.zeros((NoSTEPS,),float) # same as above
L_NH4 = np.zeros((NoSTEPS,),float) # same as above
TotN = np.zeros((NoSTEPS,),float) # same as above
# Mussels
Soma = np.zeros((NoSTEPS,),float) # makes a vector array of zeros (size: NoSTEPS rows by ONE column)
Gonad = np.zeros((NoSTEPS,),float) # same as above
B = np.zeros((NoSTEPS,),float) # Biomass - same as above
L_Temp = np.zeros((NoSTEPS,),float) # same as above
L_Salt = np.zeros((NoSTEPS,),float) # same as above
L_Oxy = np.zeros((NoSTEPS,),float) # same as above
L_Food = np.zeros((NoSTEPS,),float) # same as above
F = np.zeros((NoSTEPS,),float) # same as above
A = np.zeros((NoSTEPS,),float) # same as above
R = np.zeros((NoSTEPS,),float) # same as above
RE = np.zeros((NoSTEPS,),float) # same as above
Spawning = np.zeros((NoSTEPS,),float) # same as above
Copy-paste the section "Initializing with initial conditions" from Mussels.py to NPZD2_Mussels.py
Copy-paste the "Mussel" section just below the one for "NPZD2".
The new code should look like this:
In [ ]:
# Initializing with initial conditions
Phy[0] = InitCond['Phy']
Zoo[0] = InitCond['Zoo']
SDet[0] = InitCond['SDet']
LDet[0] = InitCond['LDet']
NH4[0] = InitCond['NH4']
NO3[0] = InitCond['NO3']
# Mussels
Soma[0] = InitCond['Soma']
Gonad[0] = InitCond['Soma']
B[0] = InitCond['Soma'] + InitCond['Gonad']
Spawning[0] = 0.
Salt = np.ones((NoSTEPS,),float) * InitCond['Salt']#Salinity
Oxy = np.ones((NoSTEPS,),float) * InitCond['Oxy'] #Oxygen
Copy-paste "Update and step" from Mussels.py to NPZD2_Mussels.py Copy-paste the "Mussel" "Update and step" just below the one for "NPZD2".
The new code should look like this:
In [ ]:
# Update and step ----------------------------------------------------
Phy[t+1] = Phy[t] + (dPhydt * dt)
Zoo[t+1] = Zoo[t] + (dZoodt * dt)
SDet[t+1] = SDet[t] + (dSDetdt * dt)
LDet[t+1] = LDet[t] + (dLDetdt * dt)
NH4[t+1] = NH4[t] + (dNH4dt * dt)
NO3[t+1] = NO3[t] + (dNO3dt * dt)
if NO3[t+1] <= 0.0001:
offset = NO3[t+1]
NH4[t+1] = NH4[t+1] + offset
NO3[t+1] = NO3[t+1] - offset
# Mussels
Soma[t+1] = Soma[t] + (dSomadt * dt)
Gonad[t+1] = Gonad[t] + (dGonaddt * dt) - Spawning[t]
B[t+1] = Soma[t+1] + Gonad[t+1]
Copy-paste MAIN MODEL LOOP from Mussels.py to NPZD2_Mussels.py Copy-paste the "Mussel" MAIN MODEL LOOP just below the one for "NPZD2"
The new code should look like this:
In [ ]:
# *****************************************************************************
# MAIN MODEL LOOP *************************************************************
for t in range(0,NoSTEPS-1):
muMax = par['mu0'] * (1.066 ** Temp[t]) # text
f_I[t] = (par['alpha']*I[t])/(np.sqrt(muMax**2+((par['alpha']**2)*(I[t]**2)))) #Eq5
L_NO3[t] = max(0,((NO3[t])/(par['kNO3']*NO3[t])) * (1/(1+(NH4[t]/par['kNH4'])))) #Eq3
L_NH4[t] = max(0,NH4[t]/(par['kNH4']+NH4[t])) # Eq 4
mu[t] =muMax * f_I[t] * (L_NO3[t] + L_NH4[t]) # Eq2
g = par['gmax'] * (Phy[t]**2/(par['kP']+(Phy[t]**2)))
n = par['nmax'] * (1 - max(0,(I[t]-par['I0'])/(par['kI']+I[t]-par['I0'])))
dPhydt = (mu[t] * Phy[t]) - \
(g * Zoo[t]) - \
(par['mP'] * Phy[t]) - \
(par['tau']*(SDet[t]+Phy[t])*Phy[t]) # Eq1
dZoodt = (g * par['beta'] * Zoo[t]) - \
(par['lBM']*Zoo[t]) - \
(par['lE']*((Phy[t]**2)/(par['kP']+(Phy[t]**2)))*par['beta']*Zoo[t]) - \
(par['mZ']*(Zoo[t]**2))#Eq10
dSDetdt = (g * (1-par['beta']) * Zoo[t]) + \
(par['mZ']*(Zoo[t]**2)) + \
(par['mP'] * Phy[t]) - \
(par['tau']*(SDet[t]+Phy[t])*SDet[t]) - \
(par['rSD']*SDet[t])
dLDetdt = (par['tau']*((SDet[t]+Phy[t])**2)) - \
(par['rLD']*LDet[t])
dNO3dt = -(muMax * f_I[t] * L_NO3[t] * Phy[t]) + \
(n * NH4[t])
dNH4dt = -(muMax * f_I[t] * L_NH4[t] * Phy[t]) - \
(n * NH4[t]) + \
(par['lBM'] * Zoo[t]) + \
(par['lE']*((Phy[t]**2)/(par['kP']+(Phy[t]**2)))*par['beta']*Zoo[t]) + \
(par['rSD']*SDet[t]) + \
(par['rLD']*LDet[t])
# MUSSELS -------------------------------------------------------------
# Calculate Temperature Limitation
L_Temp[t] = min(max(0.,1.-np.exp(-par['KTempL']*(Temp[t]-par['TempL']))), \
max(0.,1.+((1.-np.exp(par['KTempH']*Temp[t]))/(np.exp(par['KTempH']*par['TempH'])-1.))))
# Calculate Salinity Limitation
L_Salt[t] = max(0.,1.-np.exp(-par['KSaltL']*(Salt[t]-par['SaltL'])))
# Calculate Oxygen Limitation
L_Oxy[t] = max(0.,1.-np.exp(-par['KOxyL']*(Oxy[t]-par['OxyL'])))
# Calculate Oxygen Limitation
L_Food[t] = (Phy[t]+Zoo[t]+SDet[t])/(par['KFood']+Phy[t]+Zoo[t]+SDet[t])
# Calculate Filtration rate
Fmax = par['Fmax_ref']*(B[t]**(2./3.))
F[t] = Fmax * L_Temp[t] * L_Salt[t] * L_Oxy[t] * L_Food[t]
A[t] = F[t] * ((par['epsilonP']*par['AE_P']*Phy[t])+ \
(par['epsilonZ']*par['AE_Z']*Zoo[t])+ \
(par['epsilonD']*par['AE_D']*SDet[t]))
R[t] = (par['Rm']*B[t]) + (par['beta']*A[t])
RE[t] = max(0., (B[t]-par['Bpub'])/(par['KRE'] + B[t] - (2.*par['Bpub'])))
# Spawning
if Gonad[t]/B[t] < par['GT']:
Spawning[t] = 0.
dGonaddt = (A[t]-R[t]) * RE[t]
dSomadt = (A[t]-R[t]) * (1.-RE[t])
elif Gonad[t]/B[t] >= par['GT']:
Spawning[t] = Gonad[t]
dGonaddt = 0.
dSomadt = A[t]-R[t]
Copy-paste section "Pack to output" from Mussels.py to NPZD2_Mussels.py Copy-paste the "Mussel" section just below the one for "NPZD2". Note that they are duplicates, son;t copy the duplicates (just Soma, Gonad, B, and Spawning).
The new code should look like this:
In [ ]:
# Pack output into dictionary
output = {}
output['time'] = time
output['Phy'] = Phy
output['Zoo'] = Zoo
output['SDet'] = SDet
output['LDet'] = LDet
output['NH4'] = NH4
output['NO3'] = NO3
output['mu'] = mu
output['f_I'] = f_I
output['L_NO3'] = L_NO3
output['L_NH4'] = L_NH4
output['TotN'] = TotN
output['Soma'] = Soma
output['Gonad'] = Gonad
output['B'] = B
output['Spawning'] = Spawning
Update plot by adding to the plotting code for ax2:
ax2.plot(output['time']/365,output['B'],'r.')The new code should look like this:
In [ ]:
ax2.plot(output['time']/365,output['Phy'],'g-')
ax2.plot(output['time']/365,output['Zoo'],'r-')
ax2.plot(output['time']/365,output['SDet'],'k-')
ax2.plot(output['time']/365,output['LDet'],'k-.')
ax2.plot(output['time']/365,output['NH4'],'m-')
ax2.plot(output['time']/365,output['NO3'],'c-')
ax2.plot(output['time']/365,output['TotN'],'y-')
ax2.plot(output['time']/365,output['B'],'r.')
ax2.set_xlabel('Time (years)')
ax2.set_ylabel('Nitrogen (mmol N m$^{-3}$)')
ax2.set_title('Fennel et al 2006 Model')
plt.legend(['Phy','Zoo','SDet','LDet','NH4','NO3','TotN','B'])
Run the model so that you can see a plot.
Anything changed?
Is total nitrogen conserved?
Are you sure?
"B" is not accounted for in the Total Nitrogen calculation.
Add "B[t+1]" to the "Total Nitrogen calculation. The new code should look like this:
In [ ]:
# Estimate Total Nitrogen
TotN[t+1] = Phy[t+1] + Zoo[t+1] + SDet[t+1] + LDet[t+1] + NH4[t+1] + NO3[t+1] + B[t+1]
Run the model to plot again
We added Mussels to NPZD2, but they are not connected yet. "Food" is eaten by the mussels but are not taken out of the NPZD2 pool.
To connect them, we have to add extra terms to our original equations. Below, as an example, es the original Phytoplankton equation and then the new equation with an added term to account for mussel filtration.
(original from Fennel et al. 2006)
$$ \frac{\partial Phy}{\partial t} = \mu Phy - gZoo - m_P Phy - \tau (SDet -Phy)Phy$$(new) $$ \frac{\partial Phy}{\partial t} = \mu Phy - gZoo - m_P Phy - \tau (SDet -Phy)Phy - F \epsilon_P Phy$$
To actually code the "feedback", we need to add code to Ad Feaces to LDet, eat Phy/Zoo/SDet... and excrete into NH4. The new code goes after "Spawning" but before "Update and setep. Also, we need to add spawning to Zoo.
The new code should look like this:
In [ ]:
# Spawning
if Gonad[t]/B[t] < par['GT']:
Spawning[t] = 0.
dGonaddt = (A[t]-R[t]) * RE[t]
dSomadt = (A[t]-R[t]) * (1.-RE[t])
elif Gonad[t]/B[t] >= par['GT']:
Spawning[t] = Gonad[t]
dGonaddt = 0.
dSomadt = A[t]-R[t]
#Feedback to NPZD2 model
# Faeces and Pseudofaeces
Fae = F[t] * ((par['epsilonP']*(1-par['AE_P'])*Phy[t])+ \
(par['epsilonZ']*(1-par['AE_Z'])*Zoo[t])+ \
(par['epsilonD']*(1-par['AE_D'])*SDet[t]))
dLDetdt = dLDetdt + Fae
# Remove eaten Phy, Zoo and SDet from water-column
dPhydt = dPhydt-(F[t] *par['epsilonP']*Phy[t])
dZoodt = dZoodt-(F[t] *par['epsilonZ']*Zoo[t])
dSDetdt = dSDetdt -(F[t] *par['epsilonD']*SDet[t])
# Excretion into Ammonia
dNH4dt = dNH4dt + R[t]
# Update and step ----------------------------------------------------
Phy[t+1] = Phy[t] + (dPhydt * dt)
Zoo[t+1] = Zoo[t] + (dZoodt * dt) + Spawning[t]
SDet[t+1] = SDet[t] + (dSDetdt * dt)
LDet[t+1] = LDet[t] + (dLDetdt * dt)
NH4[t+1] = NH4[t] + (dNH4dt * dt)
NO3[t+1] = NO3[t] + (dNO3dt * dt)
if NO3[t+1] <= 0.0001:
offset = NO3[t+1]
NH4[t+1] = NH4[t+1] + offset
NO3[t+1] = NO3[t+1] - offset
# Mussels
Soma[t+1] = Soma[t] + (dSomadt * dt)
Gonad[t+1] = Gonad[t] + (dGonaddt * dt) - Spawning[t]
B[t+1] = Soma[t+1] + Gonad[t+1]
Run the model to plot again
Is total nitrogen conserved?
How many time mussels spawn in a summer?
Is phytoplankton eaten by the "spawned" mussel larvae? ...can you see it in the Phy series?
In [ ]: