The nonlinear diffusion equation $\rho_t-\Delta(\rho H(\rho-\rho_c))=0$ in $(0,\infty)\times \mathbb{R}^n$, $\rho(0,x)=\rho_0(x)$, where $H$ is the Heaviside function, describes the sandpile model [P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. A (3), pp. 364--374; P. Bantay and I. M. Janosi, Phys. A, 185 (1992), pp. 11--18; R. Cafiero et al., Europhys. Lett., 29 (1995), pp. 111--116] with critical state $\rho_c\in L^\infty(\mathbb{R}^n)\cap L^1(\mathbb{R}^n)$. Here, one proves that a solution $\rho=\rho(t,x)$ can be obtained as the limit of the time-stepping approximation scheme associated with the variational problem $\rho_k={\rm arg}\,\min_{\rho\in\mathcal{P}}\{\frac1h\ d^2(\rho_{k-1},\rho)+\mathbb{E}(\rho)\}$, where $d$ is the $2$-Wasserstein distance and $\mathbb{E}$ is the energy functional corresponding to the above nonlinear diffusion process. This result is on the line of that previously obtained for the linear Fokker--Planck and porous media equations by Jordan, Kinderlehrer, and Otto [SIAM J. Math. Anal., 29 (1998), pp. 1--17] and Otto [Arch. Ration. Mech. Anal., 141 (1998), pp. 63--103; Comm. Partial Differential Equations, 26 (2001), pp. 101--174].