In this article we study the numerical solution of the $L^1$-Optimal Transport Problem on 2D surfaces embedded in $R^3$, via the DMK formulation introduced in [FaccaCardinPutti:2018]. We extend from the Euclidean into the Riemannian setting the DMK model and conjecture the equivalence with the solution Monge-Kantorovich equations, a PDE-based formulation of the $L^1$-Optimal Transport Problem. We generalize the numerical method proposed in [FaccaCardinPutti:2018,FaccaDaneriCardinPutti:2020] to 2D surfaces embedded in $REAL^3$ using the Surface Finite Element Model approach to approximate the Laplace-Beltrami equation arising from the model. We test the accuracy and efficiency of the proposed numerical scheme, comparing our approximate solution with respect to an exact solution on a 2D sphere. The results show that the numerical scheme is efficient, robust, and more accurate with respect to other numerical schemes presented in the literature for the solution of ls$L^1$-Optimal Transport Problem on 2D surfaces.