## Appendix A
### Stability Analysis  
To apply the method of frozen coefficients, we make the following assumptions. We set the diffusion coefficient to its maximum value \(D=D_{\rm max}=0.65\). Furthermore, notice that the source term \(\rho c(1-\frac{c}{c_{\rm max}})\) has a maximum function value \(C_{0}=\rho c_{\rm max}/4\) with the constant \(\rho\). We note that, by freezing these coefficients, we practically linearize the model for which the standard stability analysis is applicable. For simplicity purposes, we consider the one-dimensional version of Equation (1). The extension of the same analysis to the two dimensional model is straightforward. Finally, we can ignore the constant source term \(C_{0}\), because, by a simple transformation \(\tilde{c}=c-C_{0}\), the in-homogeneous equation becomes a homogeneous one. Now, we are ready to apply the procedures of the stability analysis [45; 46] to the following simplified model with prescribed boundary conditions  
\[c_{t}\,=\,D_{\rm max}\,c_{xx}\,,\qquad c_{x}|_{\,x=0}\,=\,0,\qquad c_{x}|_{\, x=50}\,=\,0. \tag{1}\]We consider the second-order Crank-Nicolson scheme below (the analysis for the first-order scheme would be similar)  
\[\frac{c_{i}^{n+1}-c_{i}^{n}}{k}\;=\;\frac{D_{\max}}{2}\left(\frac{c_{i+1}^{n+1}- 2c_{i}^{n+1}+c_{i}^{n+1}}{h^{2}}\right)+\;\frac{D_{\max}}{2}\left(\frac{c_{i+1}^ {n}-2c_{i}^{n}+c_{i}^{n}}{h^{2}}\right) \tag{11}\]  
where \(h\;=\;\Delta x\) and \(k\;=\;\Delta t\). Letting \(r\;=\;\frac{D_{\max}\Delta t}{\Delta x^{2}}\) (11) becomes  
\[-\frac{r}{2}c_{i-1}^{n+1}+(1+r)c_{i}^{n+1}-\frac{r}{2}c_{i+1}^{n+1}\;=\;\frac{ r}{2}c_{i-1}^{n}+(1-r)c_{i}^{n}+\frac{r}{2}c_{i+1}^{n},\;\;\;i=0,\ldots,M-1. \tag{12}\]  
Our second-order boundary condition treatment at the left boundary leads to  
\[\frac{c_{0}^{n,n+1}-c_{-1}^{n,n+1}}{h}\;=\;0\Rightarrow c_{-1}^{n,n+1}\;=\;c _{0}^{n,n+1}, \tag{13}\]  
and at the right boundary we have  
\[\frac{c_{M}^{n,n+1}-c_{M-1}^{n,n+1}}{h}\;=\;0\Rightarrow c_{M}^{n,n+1}\;=\;c _{M-1}^{n,n+1}. \tag{14}\]  
Making use of (13) and (14) in (12), we have  
\[(1+\frac{r}{2})c_{0}^{n+1}-\frac{r}{2}c_{1}^{n+1}\;=\;(1-\frac{r}{2})c_{0}^{n} +\frac{r}{2}c_{1}^{n}, \tag{15}\]  
\[-\frac{r}{2}c_{M-2}^{n+1}+(1+\frac{r}{2})c_{M-1}^{n+1}\;=\;\frac{r}{2}c_{M-2} ^{n}+(1-\frac{r}{2})c_{M-1}^{n}. \tag{16}\]  
Now, considering (12), (15) and (16), the numerical scheme in matrix form becomes  
\[\underbrace{\left[ \begin{matrix}1+\frac{r}{2}&-\frac{r}{2}&0&0&\ldots&0\\ -\frac{r}{2}&(1+r)&-\frac{r}{2}&0&\ldots&0\\ &\vdots&&\\ 0&\ldots&-\frac{r}{2}&(1+r)&-\frac{r}{2}&0\\ &\ldots&0&0&-\frac{r}{2}&1+\frac{r}{2}\end{matrix} \right]}_{Q_{1}}\times \left[ \begin{matrix}c_{0}^{n+1}\\ c_{1}^{n+1}\\ \vdots\\ c_{1}^{n+1}\\ \vdots\\ c_{M-2}^{n+1}\\ c_{M-1}^{n+1}\end{matrix} \right]}_{\mathbf{c}^{n+1}}=\]  
\[\underbrace{\left[ \begin{matrix}1-\frac{r}{2}&\frac{r}{2}&0&0&\ldots&0\\ \frac{r}{2}&1-r&\frac{r}{2}&0&\ldots&0\\ &\vdots&&\\ 0&\ldots&0&\frac{r}{2}&1-r&\frac{r}{2}\\ 0&\ldots&0&0&\frac{r}{2}&1-\frac{r}{2}\end{matrix} \right]}_{\mathbf{Q}}\times \left[ \begin{matrix}c_{0}^{n}\\ c_{1}^{n}\\ \vdots\\ c_{i}^{n}\\ \vdots\\ c_{M-2}^{n}\\ c_{M-1}^{n}\end{matrix} \right]}_{\mathbf{c}^{n}},\]  
or  
\[Q_{1}\,\mathbf{c}^{n+1}\;=\;Q\,\mathbf{c}^{n}. \tag{17}\]  
Noting that \(Q\,=\,2\,I\,-\,Q_{1}\), we have \(Q_{1}\,\mathbf{c}^{n+1}\;=\;(2\,I\,-\,Q_{1})\,\mathbf{c}^{n}\) or \(\mathbf{c}^{n+1}\;=\;(2\,Q_{1}^{-1}\,-\,I)\,\mathbf{c}^{n}\) or \(\mathbf{c}^{n+1}\;=\;\widetilde{Q}\,\mathbf{c}^{n}\) with \(\widetilde{Q}\,=\,(2\,Q_{1}^{-1}\,-\,I)\). Since \(Q_{1}\) is a symmetric matrix, its eigenvalues are real, and if \(\lambda\) is an eigenvalue of \(Q_{1}\), then \(\mu\,=\,\frac{2}{\lambda}-1\) is an eigenvalue of \(\widetilde{Q}\)[45]. Moreover, \(Q_{1}^{-1}\) and therefore \(\widetilde{Q}\) are also symmetric because  
\[I=(Q_{1}\,Q_{1}^{-1})^{T}\,=\,(Q_{1}^{-1})^{T}Q_{1}^{T}\,=\,(Q_{1}^{-1})^{T}Q_{ 1}\Rightarrow(Q_{1}^{-1})^{T}\,=\,Q_{1}^{-1}. \tag{114}\]  
Here, we intend to utilize the Proposition 3.1.13 from the text [45] to prove the stability of the numerical scheme.