Derivative of the flow for ODEs on manifolds Let $\mathbf V \colon [0,T] \times \mathbb R^d \to \mathbb R^d$ (for $T>0$) be a given, bounded smooth vector field and let $\mathbf X=\mathbf X(t,x)$ be its flow, i.e. the unique solution to the initial-value problem
\begin{equation}
\begin{cases}
\frac{\partial}{\partial t} \mathbf X(t,x) = \mathbf V(t,\mathbf X(t,x)) & \text{ in } (0,T) \times \mathbb R^d \\
\mathbf X(0,x) = x \quad \text{ for all } x \in \mathbb R^d.    
\end{cases}
\end{equation}  
A well-known result in standard ODE's theory says that 
$$\tag{1}
\nabla_x \mathbf X(t,x) = \exp\bigg( \int_0^t \nabla \mathbf V(s,\mathbf X(s,x))\,ds\bigg).   
$$

Is there an analogous formula to (1) for ODEs driven by (smooth) vector fields on Riemannian manifolds? In particular, does this formula involve somehow the geometry of the Riemannian manifold? 
  A rather precise question could be: consider the $C^1$ norm of $\mathbf X$ (or even its Lipschitz constant) w.r.t. space variable $x$: does it depend on some known tensors on the manifold (e.g. curvature)?

I have gone through books in differential geometry/differential topology (e.g. Lee, Lang) and they prove that $\mathbf X$ is smooth but do not compute explicitly the derivative.
References are very much welcome. Thanks.
 A: As was pointed out by Deane Yang and Igor Khavkine in the comments, this feels like a fact that should be "looser" than Riemannian geometry.  Indeed, as I will show below, your formula makes sense in the broader setting of smooth manifolds equipped with a volume form.  I don't know a reference for this fact.
If $(M, \Omega)$ is a manifold equipped with volume form, ${\bf V}\in \operatorname{Vect}(M)$ is a vector field on $M$, and ${\bf X}_t:M\to M$ is the flow of ${\bf V}$, then:
1) One has  $J_t\in C^\infty(M)$, (i.e., a smooth  $J:\mathbb{R}\times M\to M$), with $J_0\equiv 1$, given by, 
$$J_t(x)=(\Omega_x)^{-1}\otimes \Omega_{{\bf X}_t(x)}\otimes \Lambda^n(d{\bf X}_t|_x),$$
well-defined since $\Lambda^n(d{\bf X}_t|_x):\Lambda^n( T_xM)\to \Lambda^n (T_{{\bf X}_t(x)}M)$ as Igor said.
2) One has a well-defined $\operatorname{div}({\bf V})$, namely (up to a sign, I forget whether $\pm$)
$$\operatorname{div}({\bf V})=d(\iota_{\bf V}\Omega)/\Omega.$$
So the "hopefully-an-identity" 
$$\left.\frac{d}{dt}\right\rvert_{t=0}J_t=\operatorname{div}({\bf V})$$
(for simplicity I state just the $t=0$ formula) is well-formed.  You can then prove it by choosing local co-ordinates in which $\Omega\equiv 1$.  (There is surely an intrinsic proof, too, but I can't think of one now.) 
Afterthought, 30 Oct:  I imagine the intrinsic proof consists of identifying both sides with the Lie derivative $(\mathscr{L}_{\bf V}\Omega)/\Omega$.
Update, 30 Oct, in response to questions of OP in comments: If I calculate correctly, then yes, the formulas 
\begin{align*}
\frac{dJ_t}{dt}(x)&=\operatorname{div}({\bf V})|_{{\bf X}_t(x)} J_t(x)\\
J_t(x)&=\exp\left(\int_0^t\operatorname{div}({\bf V})|_{{\bf X}_s(x)}ds\right)
\end{align*}
are also correct.
And yes, although the formula for "$\tfrac{d}{dt}\left[\det\left(d{\bf X}_t\right)\right]$" depends only on the volume form, I think one needs a Riemannian metric $g$ in order to make sense of "$\tfrac{d}{dt}\left[d{\bf X}_t\right]$" (really, the $g$-covariant derivative $D_t$ of the section $d{\bf X}_t$ of the bundle $T_x^*M\otimes TM$ along the curve ${\bf X}_t(x)$).  Probably 
$$D_t\left[d{\bf X}_t\right]|_{t=0}=\nabla {\bf V},$$
where $\nabla$ is the Levi-Civita connection of $g$.  (Please check this calculation before relying on it!)
A: I think you may be interested in this remarkable, recent paper by E. Brué and D. Semola (see in particular Theorem 3.11 which answers to your question in a much more general setting). 
