Dependence of distance function on Riemannian metric
Solution 1:
I will show that if the metrics are $C^1$ close then so are their respective distance functions. I guess the proof can be iterated and thus one should be able to show that if the metrics are $C^k$ close then also the distance functions are $C^k$ close.
As pointed out in a comment, if two metrics are $C^0$ close, then their distance functions are also $C^0$ close. Now, let $g$ be a metric with induced distance function $d$, and note that the gradient of $d$ with respect to $g$ is given by $$ \nabla_g d(x,y) = (u_1,u_2) \in T_{(x,y)} (M \times M) \simeq T_xM \oplus T_yM , $$ where $u_1 = -\gamma'(0)$, $u_2 = \gamma'(d(x,y))$, and $\gamma: [0,d(x,y)] \to M$ is the unique unit speed geodesic joining $x$ to $y$. In other words, we have $$ \nabla_g d(x,y) = (\nabla_g d^y(x),\nabla_g d^x(y)) , $$ where $d^z(\cdot)$, $z \in M$, is the "distance to $z$" function.
We argue by contradiction. Let the sequence of metrics $g_n \to g$ in $C^1$, and assume that the sequence of induced distance functions $d_n$ does not converge in $C^1$ to $d$. Since they do converge in $C^0$, this means there is a sequence $(x_n,y_n) \in M \times M$ and an $\epsilon > 0$ such that $$ |\nabla_{g_n} d_n(x_n,y_n) - \nabla_g d(x_n,y_n)|_g \geq \epsilon . \tag{1} $$ Without loss of generality assume $(x_n,y_n) \to (x,y)$, $x \neq y$, and note that $d_n(x_n,x_n) \to d(x,y)$ by $C^0$ convergence. Let $\gamma_n$ be the unique unit speed geodesics in the metrics $g_n$ joining $x_n$ to $y_n$, and let $\gamma$ be the unique geodesic joining $x$ to $y$. Write $\gamma_n'(0) = u_n$, $\gamma'(0) = u$. A simple calculation shows $$ ||u_n|_g - 1| = O(\|g_n-g\|_{C^0}) $$ and thus we can assume without loss of generality that $(x_n,u_n)$ converges to some $(x,w)$ in $TM$, where $|w|_g = 1$. By the inequality in $(1)$ we can also suppose $w \neq u$. Consider the geodesic $\alpha$ (in $g$) given by $$ \alpha(0) = y , \quad \alpha'(0) = w . $$ Since $g_n \to g$ in $C^1$, by ODE theory we have that $\gamma_n \to \alpha$ in $C^1$. Also, since $\gamma$ is the unique geodesic joining $x$ and $y$, and by our hypothesis $w = \alpha'(0) \neq \gamma'(0) = u$, we have $\alpha(s) \neq y$ for all $s \in [0, d(x,y) + \delta)$ for some $\delta > 0$. All these imply that $y_n = \gamma_n(d_n(x_n,y_n)) \to \alpha(d(x,y)) \neq y$, which is the desired contradiction.