formula we get by induction that $M^{n} \in {\mathcal{K}}_h$ for any $n=0,\dots,N$.

We can provide an alternative proof to the discrete mass conservation property $(M^n,\bar 1)_2=1$ by the following observation. Since for integrating on the torus ${\mathbb{T}}^2$ the second equation, noting $$\int_{{\mathbb{T}}^2}\Delta m-{{\rm{div}}}(mH_p(x,Du))dx=0,$$ so that $$\begin{aligned}
0&=\int_{{\mathbb{T}}^2}\partial_{(0,t]}^{\alpha}m(x,t)dx=\frac{1}{\Gamma(1-\alpha)}\int_{{\mathbb{T}}^2}\int_0^t\frac{1}{(t-s)^{\alpha}}\frac{\partial m(x,s)}{\partial s}dsdx\\
&=\frac{1}{\Gamma(1-\alpha)}\int_0^t\frac{1}{(t-s)^{\alpha}}\big{(}\int_{{\mathbb{T}}^2}\frac{\partial m(x,s)}{\partial s}dx\big{)}ds,\,\,\,\,\,\forall t\in [0,T].\end{aligned}$$ Then $$\frac{d}{dt}\int_{{\mathbb{T}}^2}m(x,t)dx=\int_{{\mathbb{T}}^2}\frac{\partial}{\partial t}m(x,t)dx=0,\,\,\,\,\forall t\in [0,T].$$ Therefore we have $$\int_{{\mathbb{T}}^2}m(x,t)dx=\int_{{\mathbb{T}}^2}m(x,0)dx=1,\,\,\,\,\forall t\in [0,T].$$ Analogously, for the discrete case we can obtain $$\sum_{i,j}(\Delta_h M^{n+1})_{i,j} +{\mathcal{B}}_{i,j}(U^{n},M^{n+1})=0$$ and $$\begin{aligned}
0=h^2\sum_{i,j}D^{\alpha}_{\Delta t}M_{i,j}^{n+1}=D^{\alpha}_{\Delta t}(M^{n+1},\bar{1})_2,\,\,\,0\leq n\leq N-1.\end{aligned}$$ By using the equivalent form of L1 scheme for forward Caputo derivative we have $$\sum_{k=0}^{n}b_k\frac{(M^{n+1-k},\bar{1})_2-(M^{n-k},\bar{1})_2}{\Delta t}=0,\,\,\,0\leq n\leq N-1.$$ It follows easily by mathematical induction that $$(M^{n},\bar{1})_2=(M^{0},\bar{1})_2=1,\,\,\,1\leq n\leq N.$$

The following lemma provides an integration by parts for , which is the discrete counterpart of fractional integration by parts formula with Caputo derivative (see [@q Lemma 2.2]) $$\int_{0}^T \partial_{(0,t]}^{\alpha} \phi(t)\kappa(t)dt + \phi(0)(I_{[t,T)}^{1-\alpha}\kappa)(0) = \int_{0}^T \phi(t) \partial_{[t,T)}^{\alpha} \kappa(t)dt + \kappa(T)(I_{(0,t]}^{1-\alpha}\phi)(T),$$ where $\phi(t), \kappa(t)\in C^1[0,T]$. \[le:by\_parts\_n\_n\_plus\_1\] Given $U=\{U^n\}_{n=0}^{N}$, $M=\{M^n\}_{n=0}^{N}$, then $$\begin{aligned}
&\sum_{n=0}^{N-1} \Big(\bar D_{\Delta t}^{\alpha} U^{n}, M^{n+1}\Big)_2 + \frac{1}{{\rho_\alpha}}\sum_{n=0}^{N-1} \bar c_n^N \left(U^N,M^{n+1}\right)_2  \\
&= \sum_{n=0}^{N-1} (D_{\Delta t}^{\alpha} M^{n+1},U^{n})_2+ \frac{1}{{\rho_\alpha}}\sum_{n=0}^{N-1} c_0^{n+1}\left(M^0,U^n\right)_2.\end{aligned}$$ For the proof, see the Appendix. We also recall the following result [@glm Lemma 2.2], which will be used for 