Stepping Through the Nitrogen Cycle, 0.1 Hours at a Time

A reef tank’s nitrogen cycle doesn’t jump from state to state—it flows. Ammonia doesn’t vanish and become nitrite instantly; it decays continuously as Nitrosomonas bacteria convert it, molecule by molecule, while new waste keeps arriving. You can write the rate as dA/dt = -k₁·A, but you can’t solve for hour 73 without knowing what happened in hour 72.9.

Runge-Kutta methods let you step through continuous change by sampling the derivative at multiple points within each timestep, then blending those slopes into a single weighted step forward. RK4—the fourth-order variant—hits the sweet spot: four evaluations per step, error proportional to h⁵.

Here’s a minimal RK4 stepping through 10 days of a cycling tank, where ammonia decays (k₁ = 0.05/hour) into nitrite, which then decays (k₂ = 0.03/hour) into nitrate:

rk4 <- function(state, deriv, dt, steps) {
  for (i in 1:steps) {
    k1 <- deriv(state)
    k2 <- deriv(state + 0.5*dt*k1)
    k3 <- deriv(state + 0.5*dt*k2)
    k4 <- deriv(state + dt*k3)
    state <- state + (dt/6)*(k1 + 2*k2 + 2*k3 + k4)
    if (i %% 24 == 0) cat(sprintf("Day %d: NH3=%.2f NO2=%.2f NO3=%.2f\n", i/24, state[1], state[2], state[3]))
  }
}
deriv <- function(s) c(-0.05*s[1], 0.05*s[1]-0.03*s[2], 0.03*s[2])
rk4(c(4.0, 0.0, 0.0), deriv, 1, 240)

Same logic in Pascal, where you’d write this for a simulation of the exact cycling schedule I’m watching in my 30L tank:

program NitrogenCycle;
type Vec3 = array[1..3] of real;
var state: Vec3; i: integer;
function Deriv(s: Vec3): Vec3;
begin Deriv[1] := -0.05*s[1]; Deriv[2] := 0.05*s[1]-0.03*s[2]; Deriv[3] := 0.03*s[2]; end;
procedure RK4Step(var s: Vec3; dt: real);
var k1,k2,k3,k4: Vec3; j: integer;
begin k1:=Deriv(s); for j:=1 to 3 do k2[j]:=Deriv([s[1]+0.5*dt*k1[1],s[2]+0.5*dt*k1[2],s[3]+0.5*dt*k1[3]])[j];
  for j:=1 to 3 do k3[j]:=Deriv([s[1]+0.5*dt*k2[1],s[2]+0.5*dt*k2[2],s[3]+0.5*dt*k2[3]])[j];
  k4:=Deriv([s[1]+dt*k3[1],s[2]+dt*k3[2],s[3]+dt*k3[3]]);
  for j:=1 to 3 do s[j]:=s[j]+(dt/6)*(k1[j]+2*k2[j]+2*k3[j]+k4[j]); end;
begin state[1]:=4.0; state[2]:=0.0; state[3]:=0.0;
  for i:=1 to 240 do begin RK4Step(state,1.0); if i mod 24=0 then
    writeln('Day ',i div 24,': NH3=',state[1]:0:2,' NO2=',state[2]:0:2,' NO3=',state[3]:0:2); end; end.

The output climbs and falls in waves: ammonia drops, nitrite spikes around day 5, then nitrate accumulates as the endpoint. By day 10 you’re at NH3=0.24, NO2=1.01, NO3=2.75—close enough to the test kit readings that I trust the rate constants I pulled from Nitrosomonas growth curves. When you can’t wait a month to see if your biofilter works, you step through the math instead.