The nonlinear diffusion equation rho(t) - Delta(rho H(rho - rho c)) = 0 in (0,infinity) x 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) is an element of L-infinity(R-n) boolean AND L-1(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) = arg min(rho is an element of P){1/h d(2)(rho(k-1),rho) + E(rho)}, where d is the 2-Wasserstein distance and 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].