The higher-order estimates for the distance function Let $M$ be a complete Riemannian manifold such that inj$(M)\geq l>0$ and $|\nabla^k\text{Rm}|\leq A_k$ for any $k \geq 0$ . For a point $p$ on $M$, we have a distance function $r(x)=d(x,p)$. For any $k \geq 0$, can we find a constant $C_k$ which only depends on $a,b,A_i,l$, such that $|\nabla^kr(x)|\leq C_k$ provided $0<a<r(x)<b<l$?

In local harmonic coordinates about $p$, using Elliptic $L^p$ estimates and the formula
$$ \Delta_g g_{ij} = -2Ric_{ij} + Q(\partial g, g),$$ one obtains uniform $L^{k+2,p}$ estimates for $g_{ij}$, any $1 < p < \infty,$ (depending upon $l$ and $A_k$), and thus uniform $C^{k+1,\alpha}$ estimates, any $0< \alpha < 1$, by the Sobolev embedding theorem. On the other hand, so long as $0 < r(x) < l$, we have that $$\Delta_g r(x) = H(x),$$ where $H(x)$ 
is the mean curvature of the level set $\{ y \textrm{ : } r(y) = x \}$. If $g_{ij}$ is uniformly bounded in $C^{k+1,\alpha},$ it follows that $H$ is uniformly bounded in $C^{k,\alpha}$ in the same local harmonic coordinate chart. Now applying a Schauder estimate
to the above equation, we get a uniform $C^{k + 2,\alpha}$ bound on $r(x)$. Since this 
is an interior estimate, the bound does depend upon $a$ and $b$ in your problem statement, as well as the constant $A_k$. Moreover, it is valid in any arbitrarily chosen harmonic coordinate system, but may not be valid in other coordinates. Looking at the definition of $\nabla^k$ in local coordinates, perhaps it won't be hard to translate this into a coordinate invariant statement.
