clear from this special case. Let ${\mathbb{T}}^2_h$ be a uniform grid on the torus with step $h$, (this supposes that $N_h = 1/{h}$ is an integer), and denote by $x_{i,j}$ a generic point in ${\mathbb{T}}^2_h$ (an anisotropic mesh with steps $h_1$ and $h_2$ is possible too and we have taken $h_1=h_2$ only for simplicity). The values $U^n_{i,j}$ and $M^n_{i,j}$ denotes the numerical approximation of the function $u$ and $m$ at $(x_{i,j},t_n)=(ih,jh, n{{\Delta t}})$, $0\le i,j\le N_h$, $n=0,\dots,N$ (assuming that $N=T/{{\Delta t}}$ is an integer). We also denote by $U^n$ and $M^n$ the grid function taking the values $U^n_{i,j}$ and $M^n_{i,j}$ at $x_{i,j}\in {\mathbb{T}}^2_h$. For a generic pair of grid functions $U$, $V$ on ${\mathbb{T}}^2_h$, we define the scalar product\
$$(U,V)_2=h^2\sum_{  i,j } U_{i,j} V_{i,j},$$ the norm\
$$\|U\|_2^2=(U,U)_2$$ and the norm $$\|U\|_\infty=\max_{0\le i,j <N_h}|U_{ij}|.$$ We consider the compact and convex set $${\mathcal{K}}_h=\{M= (M_{i,j})_{ 0\le i,j <N_h}:   (M,\overline 1)_2=1;\quad   M_{i,j}\ge 0 \}\,,$$ where $\overline 1$ denotes the $N_h^2$-tuple with all components equal to $1$. Note that ${\mathcal{K}}_h$ can be viewed as the set of the discrete probability measures on ${\mathbb{T}}^2_h$.\
We first recall the L1 numerical method for the discretization of the forward and backward Caputo time-fractional derivatives (see [@jlz; @lx]). The 