Modelling nitrate in my aquarium

In the first post, I introduced my aquarium and explained why I monitor nitrate. The quantity I am particularly interested in is how quickly nitrate accumulates between water changes, after plant uptake and other losses.

To estimate that, I need a model for how the concentration changes over time and how my water tests relate to it. I have kept the model fairly small: three quantities, with the effects of feeding, nitrification and plant uptake combined into a net accumulation rate.

The state

The model describes the aquarium using the state vector

$$x(t)= \begin{pmatrix} N(t)\\ r(t)\\ T(t) \end{pmatrix},$$

where $N$ is nitrate concentration in the tank, $r$ is its net accumulation rate between water changes, and $T$ is nitrate concentration in the tap water. Time is measured in weeks, so $N$ and $T$ have units of parts per million (ppm) and $r$ has units of ppm/week.

The water in the tank is assumed to be well mixed, with approximately constant water volume between water changes. Therefore, the rate of change of nitrate concentration is

$$\frac{\mathrm{d}N}{\mathrm{d}t}=r.$$

A positive $r$ means nitrate is accumulating. A negative $r$ means the processes removing nitrate (e.g., plant uptake) outweigh those adding it (e.g., feeding). Even if I knew $r$ exactly, it would not tell me separately how much nitrate was produced and how much the plants removed.

Tap water belongs in the state vector because it contains nitrate and is what I use to refill the aquarium. Its nitrate concentration affects what happens at every water change.

Between water changes

If $r$ were constant, predicting nitrate would be straightforward:

$$N(t+\Delta t)=N(t)+r(t)\Delta t.$$

The net balance can change as plants grow, feeding changes or the aquarium settles after an intervention. I therefore allow $r$ to drift randomly:

$$\begin{aligned} \mathrm{d}N &= r\,\mathrm{d}t,\\ \mathrm{d}r &= \sqrt{q}\,\mathrm{d}W_r. \end{aligned}$$

Here $W_r$ is a standard Wiener process and $q$ controls how much the rate is allowed to change. Over an interval $\Delta t$, the change in $r$ has mean zero and variance $q\Delta t$. A larger $q$ permits more variation in the rate.

The random forcing represents variation in the net balance that has not been modelled explicitly. For a given current state, $x(t)$, the expected nitrate concentration follows the current slope, with uncertainty from rate changes growing over the forecast.

This is an integrated Wiener-process model, also used in target tracking as a nearly constant velocity model.¹ Nitrate concentration takes the role of position, and accumulation rate the role of velocity.

For the tap water, I assume a separate random walk,

$$\mathrm{d}T=\sqrt{q_T}\,\mathrm{d}W_T,$$

with independent noise. This allows the supply concentration to vary over time, with $q_T$ controlling the size of those changes.

For this model, the equations give an exact transition over any elapsed time between water changes:

$$x(t+\Delta t)=F(\Delta t)x(t)+\eta, \qquad \eta\sim\mathcal{N}\bigl(0,Q(\Delta t)\bigr),$$

where the new increment $\eta$ is independent of the current state. The transition matrix is

$$F(\Delta t)= \begin{pmatrix} 1 & \Delta t & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix}.$$

Mathematically, $F$ maps the current state to its conditional mean after an interval $\Delta t$. The first row gives $N+r\Delta t$, while the other two rows retain the current values of $r$ and $T$. Physically, the expected nitrate concentration follows the current net accumulation rate, with no systematic change in the rate or tap-water concentration over that interval. Actual values can vary through the random increment $\eta$.

The process-covariance matrix is

$$Q(\Delta t)= \begin{pmatrix} q\Delta t^3/3 & q\Delta t^2/2 & 0\\ q\Delta t^2/2 & q\Delta t & 0\\ 0 & 0 & q_T\Delta t \end{pmatrix}.$$

This is the covariance of $\eta$: the new uncertainty accumulated during the interval. The diagonal entries give the variances of its nitrate, rate and tap-water components. The rate component has variance $q\Delta t$, while integrating those rate changes gives the $q\Delta t^3/3$ contribution to nitrate variance. The tap-water component has variance $q_T\Delta t$.

The off-diagonal terms $q\Delta t^2/2$ describe correlated changes in nitrate and its accumulation rate. A rate that wanders upwards also tends to produce extra nitrate during the same interval. The zeros involving $T$ express the assumption of independent tap-water driving noise; they do not require the uncertain tank and tap concentrations themselves to be independent.

The forecast also carries forward uncertainty already present in the current state. Using the actual elapsed time $\Delta t$ lets the model handle irregular gaps between tests.

Water changes

Water changes are recorded as interventions. Draining and refilling the tank changes its nitrate concentration, so the model needs to account for that intervention.

Suppose the water volume is $V$, the tank concentration before draining is $N^-$, and I replace a fraction $f$ with tap water with nitrate concentration $T$. Assuming the tank is well mixed, the nitrate remaining after draining is proportional to $(1-f)VN^-$, and the replacement water contributes $fVT$.

After refilling to the original volume,

$$N^+=(1-f)N^-+fT.$$

The superscripts denote the concentrations immediately before and after the water change. The change is assumed to be instantaneous relative to the intervals between water tests. It changes $N$ directly; $r$ and $T$ are left unchanged at that instant.

For example, if the tank contains 10 ppm, the tap water contains 3 ppm, and I replace 25% of the water, the model gives

$$N^+=0.75\times10+0.25\times3=8.25\ \mathrm{ppm}.$$

A partial water change moves the concentration towards that of the incoming water. The size and direction of the change depend on the difference between the two concentrations as well as on how much water I replace.

The fraction $f$ is uncertain too: I estimate it from the drop in water level.

Writing $\bar f$ for its mean and $\sigma_f^2$ for its variance, and assuming $f$ is independent of the state, gives

$$\mathbb{E}[N^+] =(1-\bar f)\mathbb{E}[N^-]+\bar f\mathbb{E}[T]$$

and

$$\operatorname{Var}(N^+) = \operatorname{Var}\bigl((1-\bar f)N^-+\bar fT\bigr) +\sigma_f^2\,\mathbb{E}\bigl[(T-N^-)^2\bigr].$$

The first term includes the existing uncertainty in tank and tap concentrations, and their covariance. The second accounts for uncertainty in the amount of water replaced. Its effect grows with the difference between tank and tap nitrate: when their concentrations are similar, replacing a little more or less water makes little difference.

These two moments are exact under the independence assumption, but they do not determine the full distribution. For a fixed $f$, the water-change rule is linear in the state, so a jointly Gaussian starting state remains Gaussian. With uncertain $f$, the fraction multiplies uncertain concentrations. Each possible fraction can give a Gaussian distribution with a different mean and variance; averaging over those possibilities produces a mixture, which need not be Gaussian.

This matters because distributions with the same mean and variance can assign different probabilities to the next colour-chart reading, or to unusually high or low nitrate concentrations. As we shall see in the next post, the filter approximates the updated joint state by a Gaussian with the same mean and covariance. Those moments are preserved, but some information about the shape is lost. How much that affects the predictions is something to test.

What the test measures

Consumer nitrate tests typically require users to match the colour of a treated water sample to a printed chart. For my particular test, the freshwater scale is labelled

$$<0.5,\ 1,\ 3,\ 5,\ 15,\ 30,\ 60,\ 120,\ 200\ \mathrm{ppm}.$$

I record the chosen category. To relate it to the state, I assume a noisy continuous assay response,

$$z^*=c+\varepsilon, \qquad \varepsilon\sim\mathcal{N}(0,\sigma^2),$$

where $c=N$ for a tank test and $c=T$ for a tap-water test. The reported swatch tells me which interval contains $z^*$.

I choose to place the boundaries at geometric midpoints between adjacent positive labels, except that I use 0.5 ppm as the boundary above the lowest category. Under that choice, the 5 ppm swatch corresponds to an assay response between

$$\sqrt{3\times5}\approx3.87 \quad\text{and}\quad \sqrt{5\times15}\approx8.66\ \mathrm{ppm}.$$

Those boundaries apply to the noisy assay response. The true concentration can lie outside that interval and still produce the same swatch.

For a category with boundaries $a$ and $b$, the probability of observing it at concentration $c$ is

$$\operatorname{P}(y\mid c) = \Phi\left(\frac{b-c}{\sigma}\right) - \Phi\left(\frac{a-c}{\sigma}\right),$$

where $\Phi$ is the standard normal cumulative distribution function. I currently use $\sigma=2$ ppm and assume it does not depend on concentration. The outermost categories extend to minus and plus infinity in the assay-response model.

I should emphasise that the boundary locations and noise level are modelling assumptions. They determine how strongly a colour comparison should affect the estimates, so they deserve checking (which I will do in a later post) along with the dynamics.

What I can infer

The model keeps track of a net rate. An inferred reduction in that rate could be consistent with increasing plant uptake, but changes in feeding or other processes could also account for it. Nitrate readings alone cannot identify the cause.

There are mathematical limitations too. The model permits negative values of $N$ and $T$, even though physical concentrations are non-negative. That matters near zero, and the uncertainty grows without bound over a long forecast. Its intended use is over the short intervals between water tests.

The remaining step is to update the distribution of $N$, $r$ and $T$ as each new tank or tap-water result arrives. I use a Gaussian assumed-density filter for that, which I will work through in the next post.

Footnotes

1. Simo Särkkä and Arno Solin (2019), Applied Stochastic Differential Equations, Cambridge University Press, Section 6.2, Example 6.3. ↩