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][1] 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.
1. Using the $C^0$-scale control, we show that $d$ and $\delta_0$ are $C^0$ close.
2. By bounding the total acceleration of $g$-geodesics in the $x$-coordinates, we show that $d$ and $\delta_0$ are $C^1$ close.
3. By bounding the total jerk of $g$-geodesics in the $x$-coordinates, we show that $d$ and $\delta_0$ are $C^2$ close.
4. 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.
[1]: https://arxiv.org/abs/2001.06155