for , which relies on the monotonicity of $g$ and the positivity of the coefficients ${\overline c}_n^k$.

To prove that a solution of is uniformly Lipschitz continuous, we need the following fractional version of the discrete Gronwall Lemma (see [@lwz Lemma 3.1]).

Suppose that the nonnegative sequences $\{{\omega}^n,g^n\}_{n\in{\mathbb{N}}}$ satisfy $D_{\Delta t}^\alpha{\omega}^1\le {\lambda}_1 {\omega}^1+{\lambda}_2 {\omega}^0+g^1$ and $$D_{\Delta t}^\alpha {\omega}^n\le {\lambda}_1{\omega}^n+{\lambda}_2{\omega}^{n-1}+{\lambda}_3{\omega}^{n-2}+g^n,\quad n\ge 2$$ where ${\lambda}_1,{\lambda}_2,{\lambda}_3$ are given nonnegative constants independent of ${{\Delta t}}$. Then, there exists $\delta>0$ such for ${{\Delta t}}<\delta$, $$\label{eq:discrete_Gronwall}
{\omega}^n\le \left({\omega}^0+\frac{(n{{\Delta t}})^\alpha}{\Gamma(1+\alpha)}\max_{1\le j \le n}g^j\right)E_\alpha(2\lambda (n{{\Delta t}})^\alpha ),\qquad 1\le n\le N,$$ where $E_\alpha(z)=\sum_{k=0}^\infty z^k/\Gamma(1+k\alpha)$ is the Mittag-Leffler function and ${\lambda}={\lambda}_1+{\lambda}_2/(2-2^{1-\alpha})+{\lambda}_3/(2\cdot2^{1-\alpha}-1^{1-\alpha}-3^{1-\alpha})$.

\[prop:lip\_U\] Let $U=(U^n)_{n=0,\dots,N}$ be the solution of . Then, for ${{\Delta t}}$ sufficiently small, there exists a positive constant $L$ such that $$\label{eq:lip_est}
    \|D_hU^n\|_\infty \le L, \qquad  n=0,\dots,N.$$

Rewrite as the stationary equation $$U_{i,j}^{n} - \rho_\alpha (\Delta_h U^{n})_{i,j} + \rho_\alpha g(x_{i,j},[D_hU^{n}]_{i,j}) - \rho_\alpha f_h[M^{n+1}]_{i,j}= \sum_{k=n+1}^{N} {\overline c}_{n}^k U_{i,j}^k$$ and let $\Psi$ be the map $U_{i,j}^n=\Psi( \sum_{k=n+1}^{N} {\overline c}_{n}^k U_{i,j}^k)$ that associates to right hand side of the previous equation the corresponding solution. For $V \in {\mathbb{R}}^{N_h}$ and $\lambda\in{\mathbb{R}}$, we have $\Psi( V+\lambda)=\Psi( V)+\lambda$ (where, on the left side of the previous identity, $\lambda$ is interpreted as the constant function 