## 3 Numerical Scheme  
This section is devoted to the construction of a robust numerical scheme to approximate a solution to (4). This section provides an outline of the numerical scheme and ommites the proof of convergence. In contrast to [5; 42] no finite-difference scheme will be discussed, but a finite-element discretization, which allows for us to employ augmented Lagrangian methods, cf. [25; 26], on the continuous level of the model. In this section, a numerical scheme for a generalized problem in the form (3) will be discussed, by adding additional reaction terms \(\gamma_{S}c_{S}\). The respective generalized least-squares functional is given through  
\[LS\big{(}(v_{S},q_{S})_{S\in\{A,B,C\}}\big{)}:=\] \[\big{\|}D_{A}\mathrm{div}q_{A}+\gamma_{A}v_{A}+k_{2}v_{C}-k_{1}v_ {A}v_{B}+f_{A}\big{\|}_{L^{2}(\Omega)}^{2}+\big{\|}q_{A}-\
abla v_{A}\big{\|}_ {L^{2}(\Omega)}^{2}\] \[\quad+\big{\|}D_{B}\mathrm{div}q_{B}+\gamma_{B}v_{B}+k_{2}v_{C}- k_{1}v_{A}v_{B}+f_{B}\big{\|}_{L^{2}(\Omega)}^{2}+\big{\|}q_{B}-\
abla v_{B} \big{\|}_{L^{2}(\Omega)}^{2}\] \[\quad+\big{\|}D_{C}\mathrm{div}q_{C}+\gamma_{C}v_{C}-k_{2}v_{C}+ k_{1}v_{A}v_{B}+f_{C}\big{\|}_{L^{2}(\Omega)}^{2}+\big{\|}q_{C}-\
abla v_{C} \big{\|}_{L^{2}(\Omega)}^{2}\] \[\quad+\sum_{S\in\{A,B,C\}}\left[\big{\|}v_{S}-g_{S,D}\big{\|}_{L^ {2}(\Gamma_{D})}^{2}+\big{\|}\
abla v_{S}\cdot\
u_{\Gamma_{N}}-g_{S,N}\big{\|} _{L^{2}(\Gamma_{D})}^{2}\right]\]  
In contrast to the usual numerical schemes, as described above, the current augmented Lagrangian methods work on the continuous level of the problem. Applying the algorithms given in [26] and in [25] to (4) yields the following Algorithm 1:  
```
Input :\(j=0\), \(\big{(}x^{(0)},\mu^{(0)}\in X\times L^{2}(\Omega)^{\times 3}\big{)}\), \(0<\lambda^{\max}\), \(0<\alpha_{0}\), \(\tau\in]0,1[,1<\gamma,\varepsilon>0\). Output : Approximate solution \(\big{(}(c_{S,LS}^{(j+1)},P_{S,LS}^{(j+1)})_{S\in\{A,B,C\}}\big{)}\in X\) to (4) and approximation \(\lambda^{(j)}\) of the Lagrangian multiplier. Follow the steps S1: Define \(\lambda^{(j)}:=\min\big{(}\mu^{(j)},\lambda^{\max}\big{)}\) and approximate a solution to \[\big{(}(c_{S,LS}^{(j+1)},P_{S,LS}^{(j+1)})_{S\in\{A,B,C\}}\big{)}= \operatorname*{argmin}_{(v_{S},q_{S})_{S\in\{A,B,C\}}\in X}LS_{\alpha_{j}} \bigg{(}(v_{S},q_{S})_{S\in\{A,B,C\}},\frac{\lambda^{(j)}}{\alpha_{j}}\bigg{)}\] S2: Define \[\mu^{(j+1)}:=\big{(}\lambda^{(j)}-\alpha_{j}[c_{A}^{(j+1)},c_{B}^{(j+1)},c_{C }^{(j+1)}]\big{)}_{+}.\] If \(j=0\) or if \[\bigg{\|}\min\bigg{[}[c_{A}^{(j+1)},c_{B}^{(j+1)},c_{C}^{(j+1)}]_{+}\frac{ \lambda^{(j)}}{\alpha_{j}}\bigg{]}\bigg{\|}_{L^{2}(\Omega)}\leq\tau\bigg{\|} \min\bigg{[}[c_{A}^{(j)},c_{B}^{(j)},c_{C}^{(j)}],\frac{\lambda^{(j-1)}}{\alpha _{j-1}}\bigg{]}\bigg{\|}_{L^{2}(\Omega)}\] (6) then set \(\alpha_{j+1}:=\alpha_{j}\), otherwise set \(\alpha_{j+1}=\gamma\alpha_{j}\). S3: If \[\big{\|}[c_{A}^{(j+1)},c_{B}^{(j+1)},c_{C}^{(j+1)}]\cdot\mu^{(j+1)}\big{|}\leq\varepsilon\] and \[\big{\|}[c_{A}^{(j+1)},c_{B}^{(j+1)},c_{C}^{(j+1)}]-\big{(}[c_{A}^{(j+1)},c_{B}^ {(j+1)},c_{C}^{(j+1)}]\big{)}_{+}\big{\|}_{L^{2}(\Omega)}\leq\varepsilon\] then break the algorithm. S4: Set \(j=j+1\) and go to S1. ```  
**Algorithm 1**augmented Lagrangian algorithm For \(f\in L^{2}(\Omega)\), the map \((\bullet)_{+}:L^{2}(\Omega)\to L^{2}(\Omega)\) is defined as \((f)_{+}:=\max(0,f)\), almost everywhere.