algal-growth_model_P.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Tid (dage i året)
t = np.linspace(0, 365, 366)

# Parametre
mu_max = 0.5    # dag^-1
K_s = 0.02      # mg/L
m = 0.1         # dag^-1
Y = 1           # mg alger/mg P
k_w = 0.1       # m^-1
k_a = 0.05      # m^2/mg
z = 5.0         # m (halvdelen af vanddybden)
K_I = 50        # µE/m^2·s
tau = 100.0     # dage (opholdstid)
P_in = 0.05     # mg/L (fosforkoncentration i tilført vand)
A_in = 0.0      # mg/L (algekoncentration i tilført vand)

# Initialbetingelser
A0 = 0.1        # mg/L (initial algekoncentration)
P_uorg0 = 1.0   # mg/L (initial uorganisk fosfor)
y0 = [A0, P_uorg0]

# Incident lysintensitet ved overfladen som funktion af tid
def I_0(t):
    # Lysintensiteten topper midt på året (dag 182)
    return 200 + 150 * np.sin(2 * np.pi * (t - 91.25) / 365)

# Differentialligninger
def model(y, t):
    A, P_uorg = y
    # Næringsstofbegrænsning
    f_P = P_uorg / (K_s + P_uorg)
    # Lysdæmpningskoefficient
    k_d = k_w + k_a * A
    # Effektiv lysintensitet
    I_eff = I_0(t) * np.exp(-k_d * z)
    # Lysbegrænsning
    f_I = I_eff / (I_eff + K_I)
    # Specifik væksthastighed
    mu = mu_max * f_P * f_I
    # Algekoncentrationens ændring
    dAdt = mu * A - m * A - (A / tau)
    # Uorganisk fosfors ændring
    dPdt = -(mu / Y) * A + (1 / tau) * (P_in - P_uorg)
    return [dAdt, dPdt]

# Løsning af differentialligningerne
sol = odeint(model, y0, t)
A = sol[:, 0]
P_uorg = sol[:, 1]

# Beregning af effektiv lysintensitet til plot
I_eff_values = [I_0(time) * np.exp(-(k_w + k_a * A_val) * z) for time, A_val in zip(t, A)]

# Plotning af resultaterne
plt.figure(figsize=(12, 8))

# Algekoncentration
plt.subplot(3, 1, 1)
plt.plot(t, A, 'g', label='Algekoncentration (A)')
plt.ylabel('Koncentration (mg/L)')
plt.title('Årsvariation af algekoncentration, uorganisk fosfor og lysintensitet')
plt.legend()

# Uorganisk fosfor
plt.subplot(3, 1, 2)
plt.plot(t, P_uorg, 'b', label='Uorganisk fosfor (P_uorg)')
plt.ylabel('Koncentration (mg/L)')
plt.legend()

# Effektiv lysintensitet
plt.subplot(3, 1, 3)
plt.plot(t, I_eff_values, 'orange', label='Effektiv lysintensitet (I_eff)')
plt.xlabel('Tid (dage)')
plt.ylabel('Lysintensitet (µE/m²·s)')
plt.legend()

plt.tight_layout()
plt.show()