The jump-growth equation
As was observed in (Datta, Delius, and Law 2010), the variability in prey size leads to a diffusion term in the PDE for the abundance density \(N(w,t)\). This term was neglected in mizer in the past but can now be included. It becomes important when matching mizer abundances to observed abundances, because it will account for fish that are larger than average. The PDE including the diffusion term is \[ \frac{\partial N}{\partial t} = \frac12 \frac{\partial^2}{\partial w^2}(d N) - \frac{\partial}{\partial w}(g N)-\mu N \tag{1}\] where \(g\) is the growth rate, \(\mu\) is the death rate and \(d\) is a new diffusion rate.
The diffusion rate has an expression that is similar to that of the rate \(A(w)\) at which prey biomass is assimilated and thus available for metabolism, growth and reproduction. Recall that this rate is given by \[ \begin{split} A(w)=&(1-f(w))\,\gamma(w)\int N_c(w_p)\,\phi(w,w_p)\,\alpha\, w_p\,dw_p \end{split} \tag{2}\] where \(N_c(w)\) is the abundance density of prey, \(\phi(w,w_p)\) is the predation kernel, \(\gamma(w)\) is the search volume, \(f(w)\) is the feeding level and \(\alpha\) is the assimilation efficiency.
The factor \(\alpha\,w_p\) in the integral in Eq. 2 is the increase in assimilated biomass resulting from the ingestion of an individual prey of weight \(w_p\). In the expression for the diffusion rate \(D(w)\) this factor is squared: \[ D(w) = (1-f(w))\,\gamma(w) \int N_c(w_p)\phi(w/w_p)(\alpha\,w_p)^2\,dw_p. \tag{3}\]
This diffusion rate captures only the randomness in the size of encountered prey. There will be many other sources of stochasticity that we do not model explicitly in mizer. They should be added as external diffusion with ext_diffusion().
The increase \(\alpha\,w_p\) in predator mass is typically only a small proportion of the predator mass because the preferred prey are typically much smaller than the predator and also \(\alpha\) is smaller than \(1\). Because this factor is squared in the expression for \(D(w)\) it might be expected that the term in the PDE involving \(D(w)\) can be safely neglected. However it is worth testing this intuition. We will do this now by first determining the diffusion rate in a model with allometric encounter and mortality rates. We will then use that to determine its effect on the slope of the juvenile spectrum. Finally we will look at the numerical solution for the steady state.
Example calculation of diffusion rate
Let us assume that the prey abundance is given by a power law: \(N_c(w)=N_0w^{-\lambda}\) and that the predation kernel is \[ \phi(w/w_p) = \exp\left(-\frac{\log(w/w_p/\beta)^2}{2\sigma^2}\right). \tag{4}\]
and hence the integral in the expression Eq. 2 for the assimilation rate becomes \[ \begin{split} I_A&:=\int N_c(w_p)\phi(w/w_p)\alpha\,w_p\,dw_p\\ &=\alpha\int w_p^{1-\lambda}\exp\left(-\frac{\log(w/w_p/\beta)^2}{2\sigma^2}\right)dw_p. \end{split} \tag{5}\] This integral can be evaluated most easily by changing integration variable to \(x=\log(w_p/w_0)\)for an arbitrary reference weight \(w_0\) and then recognising the resulting integral \[ I_A=\alpha\, w_0^{2-\lambda}\int e^{(2-\lambda)x}\exp\left(-\frac{(x-\log(w/w_0)+\log(\beta))^2}{2\sigma^2}\right)dx \tag{6}\] as a Gaussian integral. Using the general result that \[ \int e^{ax}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)dx = \sqrt{\frac{2\pi}{\sigma^2}}\exp\left(a\mu+\frac{a^2\sigma^2}{2}\right) \tag{7}\] with \(a = 2-\lambda\) and \(\mu = \log(w/w_0)-\log(\beta)\) we find that \[ \begin{split} I_A&= \alpha\, w_0^{2-\lambda}\sqrt{\frac{2\pi}{\sigma^2}}\exp\left((2-\lambda)(\log(w/w_0)-\log(\beta))+\frac{(2-\lambda)^2\sigma^2}{2}\right)\\ &=\alpha w^{2-\lambda}\sqrt{\frac{2\pi}{\sigma^2}}\beta^{\lambda - 2}\exp\left(\frac{(2-\lambda)^2\sigma^2}{2}\right) \end{split} \tag{8}\]
Assuming an allometric search volume \(\gamma(w)=\gamma w^q\) with an exponent of \(q = n - 2 + \lambda\) we obtain \[ A(w) = (1-f(w))\,\gamma \alpha \sqrt{\frac{2\pi}{\sigma^2}}\beta^{\lambda - 2}\exp\left(\frac{(2-\lambda)^2\sigma^2}{2}\right) w^{n}. \tag{9}\]
We can evaluate the diffusion rate \(D(w)\) in the same manner. The extra factor of \(w_p\) changes \(a\) from \(2-\lambda\) to \(3-\lambda\) and we obtain \[ D(w) = (1-f(w))\gamma \alpha^2 \sqrt{\frac{2\pi}{\sigma^2}}\beta^{\lambda - 3}\exp\left(\frac{(3-\lambda)^2\sigma^2}{2}\right) w^{n+1}. \tag{10}\] Comparing this to the expression Eq. 9 for \(A(w)\) we find that \[ D(w) = A(w)w\frac{\alpha}{\beta}\exp\left(\frac{(5-2\lambda)\sigma^2}{2}\right). \tag{11}\]
To get a feel for the typical magnitude of the factor let’s consider concrete values \[ f=0.6, \qquad\alpha = 0.8, \qquad\beta = 10, \qquad\sigma = 2, \qquad\lambda = 2. \tag{12}\]
Then \[ D(w) \approx 0.59 A(w)w. \tag{13}\]
However we see that the value of \(D(w)/A(w)\) is strongly influenced by the location and width of the feeding kernel. It decreases with increasing \(\beta\) and decreasing \(\sigma\) If we choose \(\beta = 100\) and \(\sigma = 1\) then \(D(w)\approx 0.01 A(w)w\). This relative decrease in the diffusion rate with increasing \(\beta\) and decreasing \(\sigma\) may explains the result from (Datta et al. 2010) about how the stability of the system depends on these parameters.
Effect on juvenile slope
In this section we will calculate the effect of the diffusion on the slope of the juvenile spectrum in the steady state. You can skip this section if you are not interested in the details of the calculation. The outcome is that the change in the juvenile slope is small. However, the change in the size spectrum of the adults is much larger, as we will see in ?@sec-diffusion-numerical.
Without diffusion we find the juvenile spectrum by solving the steady-state equation \[ \frac{\partial}{\partial w}(g(w) N(w)) =-\mu(w) N(w). \tag{14}\]
This has the solution \[ N(w) = \frac{g(w_0)}{g(w)}N(w_0)\exp\left(-\int_{w_0}^w\frac{\mu(w')}{g(w')}dw'\right). \tag{15}\]
With allometric growth and death rates \(g(w)=g_0w^n\) and \(\mu(w)=\mu_0w^{n-1}\) this gives \[ N(w)=\left(\frac{w}{w_0}\right)^{-n} N(w_0)\exp\left(-\frac{\mu_0}{g_0}\int_{w_0}^w\frac{1}{w'}dw'\right)=N(w_0)\left(\frac{w}{w_0}\right)^{-\mu_0/g_0-n}. \tag{16}\]
Thus the juvenile steady state abundance density is given by a power law with exponent \(-\mu_0/g_0-n\).
In the presence of diffusion the steady state equation becomes the second-order ODE \[ \frac12 \frac{\partial^2}{\partial w^2}(D N) - \frac{\partial}{\partial w}(g N)-\mu N=0. \tag{17}\]
We have seen in the previous section that with \(g(w)=g_0w^n\) the diffusion rate is also a power law with one extra factor of \(w\): \(D(w)=D_0w^{n+1}\). This makes the equation Eq. 17 scale invariant and hence we again expect the solution to be a power law, So we make the Ansatz \(N(w)=N(w_0)(w/w_0)^a\) with the exponent \(a\) to be determined. Substituting this Ansatz into Eq. 17 gives \[ \frac12 D_0N_0(n+1+a)(n+a)w^{n+a-1} -g_0N_0(n+a)w^{n+a-1}-\mu_0 N_0w^{n-1+a} = 0 \tag{18}\] which requires that \[ \frac12 D_0 (n+a)^2+\left(\frac12 D_0 -g_0\right)(n+a)-\mu_0=0. \tag{19}\] This is a quadratic equation for \(n+a\): \[ \frac12 D_0 (n+a)^2+\left(\frac12 D_0 -g_0\right)(n+a)-\mu_0=0. \tag{20}\] This has two solutions \[ (n+a_\pm) = \frac1{D_0}\left(g_0-D_0/2\pm\sqrt{(g_0-D_0/2)^2+2D_0\mu_0}\right) \tag{21}\] where \(a_+\) is the solution with the + sign and \(a_-\) is the solution with the - sign.
We are only interested in the solution that goes to \(a=-\mu_0/g_0-n\) when \(D_0\to 0\). This means we are interested in the solution with the - sign. To check that indeed the solution \(a_{-}\) satisfies \(\lim_{a_-\to 0}= -\mu_0/g_0-n\) it is helpful to expand the square root term around \(D_0=0\): \[ \sqrt{(g_0-D_0/2)^2+2D_0\mu_0}=g_0+\frac{-g_0+2\mu_0}{2g_0}D_0+\frac{g_0^2-(-g_0+2\mu_0)^2}{8g_0^3}D_0^2+\dots \tag{22}\] So we find \[ \begin{split} a &= -n+\frac1{D_0}\left(g_0-D_0/2-\sqrt{(g_0-D_0/2)^2+2D_0\mu_0}\right)\\ &\approx -n-\frac{\mu_0}{g_0} - \frac18\left(1-\left(-1+2\frac{\mu_0}{g_0}\right)^2\right) D_0+\dots\\ &=-n-\frac{\mu_0}{g_0} +\frac12\left(\frac{\mu_0}{g_0}\left(1-\frac{\mu_0}{g_0}\right)\right)\frac{D_0}{g_0}+\dots \end{split} \tag{23}\]
We see that when \(\mu_0=g_0\) then the first correction term to the juvenile slope vanishes. The largest increase in slope (i.e., the least negative slope) is achieved when \(\mu_0/g_0=1/2\). In that case the slope without diffusion is \(-1.25\) (assuming \(n=0.75\)). The correction term then is \(D_0/(8g_0)\). We had already seen in the previous section that \(D_0/g_0\) is typically very small, so the change in slope is also very small.
Using the example value of \(D_0/g_0\approx0.59\) and \(n=0.75\) we get a slope correction of approximately \(0.074\) from \(-1.25\) to \(-1.176\).
Using the value \(D_0/g_0\approx0.01\) we get a smaller slope correction of approximately \(0.002\) from \(-1.25\) to \(-1.248\).
When \(\mu_0/g_0>1\) then the diffusion correction is negative, i.e., the juvenile slope becomes more negative due to diffusion, meaning fewer large fish.