You function $d_p$ is a Lipschitz function (w.r.t. to the Riemannian distance) with the Lipschitz constant $1$ and it is possibly, for every $\varepsilon>0$, to $\varepsilon$-approximate it by a smooth function with the length of the gradient $\le 1 + \varepsilon$.
The proof of this statement, even if we replace 'Riemannian' by 'Finslerian', can be found in the appendix to my paper arXiv:1112.5060; http://arxiv.org/abs/1112.5060; the proof seems to be wellknown for experts so I did not pretend that the proof is new but tried to explain it in all details. In the Riemannian category, the statement was proved in D. Azagra, J. Ferrera, F. Lopez-Mesas, Y. Rangel, Smooth approximation of Lipschitz functions on Riemannian manifolds, J. Math. Anal. Appl. 326(2007) 1370–1378. I should confess that I found it complicated to read this paper but if you need it as a citation and not in order to understand the proof it could be an option