I hope this question is still of interest. It turned out that some coauthors and I needed this fact for a recent paper. As such, I've gone through and filled in the details (and fixed some errors) from my original answer. The paper is here and the relevant section is Lemma 18 in the second appendix, which starts on page 24.
I didn't try to get explicit estimates on $\epsilon(r,Q)$, but in theory it can be done by bounding each of the constants. Sorry in advance for the length of the argument; I wasn't able to find a cleaner approach which gets around inducting on the regularity of the distance function.
To begin, we fix a point $p$, which will be the center of the ball throughout. Consider two different Riemannian metrics on $U$. The first, (denoted $g$) is the Riemannian metric of interest. The second, $g_0$, is the flat metric satisfying $(g_0){}_{ij}(x) = g_{ij}(p)$ for all $ x \in U$ (i.e. we consider the metric at $p$ and do not change the components throughout $U$). We then consider two separate distance functions. The first, which we denote $d$, is the distance from $p$ in the $g$ metric. The second, which we denote $\delta_0$, is the distance from $p$ in the $g_0$ metric.
As a broad overview, we can break up this argument into 4 steps.
- Using the $C^0$-scale control, we show that $d$ and $\delta_0$ are $C^0$ close.
- By bounding the total acceleration of $g$-geodesics in the $x$-coordinates, we show that $d$ and $\delta_0$ are $C^1$ close.
- By bounding the total jerk of $g$-geodesics in the $x$-coordinates, we show that $d$ and $\delta_0$ are $C^2$ close.
- Using the implicit function theorem, we show that the level sets of $d$ are $C^2$ close to level sets of $\delta_0$ (which are ellipsoids with bounded eccentricity). This then implies that the level sets of $d$ are convex, at least at small enough scales.
Step 1. The first step is to show that $d$ and $\delta_0$ are $C^0$ close. To do this, consider the $g$-geodesics in the $g_0$ metric (and vice versa), and use the $C^0$-scale control. This shows that the distance functions control each other, or more precisely that
$$ Q^{-1} \delta_0(p,q) \leq d(p,q) \leq Q \delta_0(p,q) $$
Step 2. We now want to show that the distance functions are close in $C^1$. To do so, consider a point $q$ which is close to $p$ (either in terms of $d$ or $\delta_0$). The geodesic from $p$ to $q$ in the $g_0$ metric is a straight line whose length (in the sense of $g_0$) is $\delta_0(q)$. Furthermore, the gradient of $\delta_0$ at $q$ is of unit norm (again in the sense of $g_0$) and points in the same direction as the line segment from $p$ to $q$.
On the other hand, the unit speed geodesic for the metric $g$ from $p$ to $q$ satisfy the equations
\begin{equation} \label{geodesicequations}
\frac{d^2 \gamma^i}{ds^2} + \Gamma^{i}_{jk}\frac{d \gamma^j}{ds}\frac{d\gamma^k}{ds} = 0
\end{equation}
where
\begin{equation} \label{Christoffel}\Gamma^i_{jk} = \frac{1}{2} g^{i\ell} \left( \frac{\partial g_{j \ell}}{\partial x^k} + \frac{\partial g_{k \ell}}{\partial x^j} - \frac{\partial g_{jk}}{\partial x^\ell} \right).
\end{equation}
However, from the scale-$C^1$-control, we can estimate that
\begin{eqnarray*}
\left \vert \frac{d^2 \gamma^i}{ds^2} \right \vert &\leq & \sum_{j,k} \frac{Q}{2} \frac{3(Q-1)}{r} \left \vert \frac{d \gamma^j}{ds} \right \vert \left \vert \frac{d \gamma^k}{ds} \right \vert
\end{eqnarray*}
along the entire geodesic. In order to make the estimates readable, we will absorb any constants involving $n,r$ and $Q$ using the notation $ f_1 \lesssim f_2$ whenever $f_1 \leq C f_2$ for some constant $ C$ depending only on $n,Q,$ and $r$. We will also use the notation
$$ \left. \left \| \frac{d^k \gamma}{ds^k} \right \|_{L^1} \right \vert_{s=\tau} = \sum_{i=1}^n \left \vert \frac{d^k \gamma^i (\tau)}{ds^k} \right \vert $$
and the corresponding notation for the $L^2$ norm as well.
In this notation, we find that
$$ \left \vert \frac{d^2 \gamma^i}{ds^2} \right \vert \lesssim \left \vert \frac{d \gamma^j}{ds} \right \vert \left \vert \frac{d \gamma^k}{ds} \right \vert. $$
Summing over the $i$ index, this implies that
\begin{eqnarray*}
\left. \left \| \frac{d^2 \gamma}{ds^2} \right \|_{L^1} \right \vert_{s=\tau} &\lesssim & \left. \left \| \frac{d \gamma}{ds} \right \|_{L^2}^2 \right \vert_{s=\tau} \\
& \lesssim & \left \| 1 + \int_0^\tau \left. \left \| \frac{d^2 \gamma}{ds^2} \right \|_{L^1} \right \vert_{s=t} \,dt \right \|_{L^2}^2 \\
& \leq & 2 + 2\left(\int_0^\tau \left. \left \| \frac{d^2 \gamma}{ds^2} \right \|_{L^1} \right \vert_{s=t} \,dt \right )^2
\end{eqnarray*}
We define $F(\tau) = \int_0^\tau \left. \left \| \frac{d^2 \gamma}{ds^2} \right \|_{L^1} \right \vert_{s=t} \,dt$. In other words $F(\tau)$ is the total acceleration of a $g$ geodesic in coordinates. When written in terms of $F$, the above estimate shows that
\begin{eqnarray*}
\frac{d F}{d \tau } &\lesssim & 1 + F^2
\end{eqnarray*}
Dividing both sides by $1 + F^2$ and integrating, we find that
$$\arctan{F(\tau)} \lesssim \tau. $$
For $\tau$ small, this provides an upper bound for $F$. However, since $\gamma$ is a unit speed geodesic, this shows that for small $d$ (equivalently $\delta_0$), the acceleration (in coordinates) of $\gamma$ is very small. As a result, $\gamma$ is $C^1$-close to a line segment from $p$ to $q$ (in coordinates). Therefore, the gradient of $d$ at $q$ is close to the gradient of $\delta_0$ at $q$, which implies that the functions are $C^1$-close as well.
Step 3. We now want to show that the distance functions are $C^2$ close as well. To do so, we start by taking one derivative of the geodesic equations and use the $C^2$ scale control to bound the $C^1$ norm of the Christoffel symbols. When we do so, we get an estimate of the form
$$\left \| \frac{ d^3 \gamma}{ d s^3 } \right \|_{L^1} \lesssim \left \| \frac{ d^2 \gamma}{ d s^2 } \right \|_{L^1} \left \| \frac{ d \gamma}{ d s } \right \|_{L^1}^2 +\left \| \frac{ d \gamma}{ d s } \right \|_{L^1}^3.$$
From here, we can basically repeat Step 2 verbatim. There are other ways of controlling this term, but it's probably simplest to do it using a method we've already done. Using Young's inequality, this shows that
\begin{eqnarray*}
\left \| \frac{ d^3 \gamma}{ d s^3 } \right \|_{L^1} & \lesssim & \left \| \frac{ d^2 \gamma}{ d s^2 } \right \|_{L^1}^2 + \left \| \frac{ d \gamma}{ d s } \right \|_{L^1}^4 +\left \| \frac{ d \gamma}{ d s } \right \|_{L^1}^3 \\
\end{eqnarray*}
Using the bound on $F$ from the $C^1$ estimate, for small time we can control $\left \| \frac{ d \gamma}{ d s } \right \|$ by a constant, which implies
\begin{eqnarray*} \left \| \dfrac{ d^3 \gamma}{ d s^3 } \right \|_{L^1} & \lesssim \left \| \dfrac{ d^2 \gamma}{ d s^2 } \right \|_{L^1}^2 + 1.
\end{eqnarray*}
Using the estimate
$$ \left. \left \| \frac{ d^2 \gamma}{ d s^2 } \right \| \right \vert_{s = \tau} \leq \int_0^\tau \left . \left \| \frac{ d^3 \gamma}{ d s^3 } \right \| \right \vert_{s = t} dt, $$
we can rewrite the preceding line to obtain
\begin{eqnarray*} \left. \left \| \frac{ d^3 \gamma}{ d s^3 } \right \| \right \vert_{s = \tau} & \lesssim \left( \int_0^\tau \left . \left \| \dfrac{ d^3 \gamma}{ d s^3 } \right \| \right \vert_{s = t} dt \right) ^2 + 1.
\end{eqnarray*}
We again integrate out the differential inequality to show that $$ \int_0^\tau \left . \left \| \dfrac{ d^3 \gamma}{ d s^3 } \right \| \right \vert_{s = t} dt \lesssim \tan(\tau). $$ For sufficiently small times $\tau$, this provides a small bound on the total jerk of $\gamma$, which also bounds the point-wise acceleration of $\gamma$. As a result, $\gamma$ is $C^2$-close to a line segment from $p$ to $q$. Therefore, the gradient of $d$ is $C^1$-close to the gradient of $\delta_0$ at q, which implies that the two distance functions are $C^2$-close.
Step 4. Finally, we show that the ball $B_{g_0}(p,\epsilon)$ is convex in the $x$-coordinates for $\epsilon$ small. Consider a point $q$ with $d(p,q)$ small and the hyperplane (in $x$-coordinates) $V$ through $q$ which is perpendicular to the line segment from $p$ to $q$. Near $q$, both the functions $d$ and $\delta_0$ have gradients which are transverse to this hyperplane. As such, by the implicit function theorem we can locally find two functions $\ell_1, \ell_2:V \to \mathbb{R}$ so that $d(v,\ell_1(v)) = d(q)$ and $d(v,\ell_2(v)) = \delta_0(q)$. In other words, we use the implicit function theorem to express the level sets of $d$ and $\delta_0$ as graphs in a small neighborhood of $q$. Furthermore, we can write the derivatives of $\ell_1$ and $\ell_2$ in terms of the derivatives of $d(q)$ and $ \delta_0(q)$, respectively.
Since $d(q)$ and $ \delta_0(q)$ are $C^2$-close, it follows that $\ell_1$ and $\ell_2$ are also $C^2$-close. However, the level sets of $\delta_0$ are ellipsoids with bounded eccentricity, so it follows that $\ell_2$ is (uniformly) strongly convex.
Since $\ell_1$ is $C^2$-close to the graph of an ellipsoid, for small enough $d$ the Hessian of $\ell_1$ must be non-negative definite. This implies that $\ell_1$ is a convex function and so secants of $\ell_1$ lie entirely above the graph of $\ell_1$. In other words, if we draw a segment between two points in $\ell_1$, the intermediate points are closer to $p$ than the endpoints. Put more simply, the ball $B_g(p,\epsilon)$ is convex (at least for $\epsilon$ small enough).
Hopefully this makes sense. In my original answer, I thought it would be possible to do a similar argument with scale-$C^1$-control, but this turns out not to be the case. It might be possible to modify the scale-$C^2$-control to allow for less smooth level sets of $d$, but the assumption needs to be strong enough to ensure that the level sets are convex. At present, I'm not sure of a good way to relax scale-$C^2$-control, and this answer is more than long enough already, so I'll leave that alone.
