is equivalent to $$\label{eq:FP_0} 
    M^{n+1}_{i,j} -\rho_\alpha(\Delta_h M^{n+1})_{i,j} - \rho_\alpha{\mathcal{B}}_{i,j}(U^{n},M^{n+1})  =  \sum_{k=0}^{n} c_k^{n+1} M_{i,j}^{k}.$$ For fixed $n$, is a linear problem in the unknown $M^{n+1}\in ({\mathbb{R}}^{N_h^2})$, which can be rewritten as $$\label{eq:FP_1}
    (I+\rho_\alpha A)M^{n+1}=\sum_{k=0}^{n} c_k^{n+1} M^{k}.$$ where $A$ is the linear operator given by $A\cdot =-  \Delta_h\cdot - {\mathcal{B}}_{i,j}(U^{n},\cdot)$. We observe that, due to the monotonicity assumption (G1), the matrix corresponding to $A$ has positive diagonal entries and non positive off-diagonal entries. Hence, for ${{\Delta t}}$ sufficiently small, the matrix corresponding to the linear operator $I+\rho_\alpha A$ in is a $M$-matrix and therefore non singular. It follows that, for any $n=0,\dots,N-1$, there exists a unique solution to . Moreover, since $M^0\ge 0$ and the inverse of a $M$-matrix is non negative, we easily get by induction that $M^n\ge 0$ for $n=0,\dots,N-1$.\
To show that $(M^n,\bar 1)_2=1$ for any $n=0,\dots,N$, we first recall the following identity (see [@acd eq. (50)]) $$\label{eq:bypart_space}
\begin{aligned}
(AW,Z)_2&=\sum_{i,j}(D^+_1W)_{i,j}(D^+_1Z)_{i,j}+\sum_{i,j}(D^+_2W)_{i,j}(D^+_2Z)_{i,j}\\
&+\sum_{i,j}W_{i,j}[D_h Z]_{i,j}\cdot D_q g(x_{i,j},[D_hU]_{i,j}).
\end{aligned}$$ From the previous identity, it follows that $(AW,\bar 1)_2=0$ for any grid function $W$. Hence, taking the inner product of with the function $\bar 1$, we get that $$(M^{n+1},\bar 1)_2=\sum_{k=0}^{n } c_k^{n+1} (M^{k},\bar 1)_2.$$ Recalling that $M^0\in {\mathcal{K}}_h$ and Lemma \[lem:c\_properties\].(iii), by the previous 