Geometric derivation of the Einstein’s field equation from the Hilbert action. It is well-known that the equation for stationary solutions of the Einstein-Hilbert functional is given by the Einstein field equation (for a statement, see previous question). The standard derivation of this is through Koszul's formulae either in coordinates (for example wikipedia), or in abstract index notation (for example, in Wald's General Relativty), or in coordinate-free notation (for example, as pointed out by Thomas Richard in Besse - Einstein manifold). This approach is mainly algebraic by using the definition in terms of Koszul's formulae and then calculus in various notations. Essentially the derivation is a direct calculation without the need to even mention the manifold. 
I am wondering if there is a way to derive/interprete the statement refered to at the beginning using an alternative method which is more geometric, ie. using parallel transport or alike. The criterion for "geometric" being (a) a direct reference to the manifold is necessary; or (b) a picture, at least in principle a mental picture, can be drawn to at least carry the main idea of the derivation (of course, pictures of formulae don't count).
 A: Comment: the following is a somewhat convoluted way of deriving the
Euler-Lagrange equation using Clairaut's theorem for the volume functional and
some standard, albeit not simpler, variation formulas (all is $C^{\infty}$ in
the following). Let $g_{0}$ be a Riemannian metric and let $v$ be a symmetric
$2$-tensor. Let $g_{0,s}$ satisfy $\frac{\partial}{\partial s}|_{s=0}
g_{0,s}=v$ and $g_{0,0}=g_{0}$. With $g_{0,s}$ as initial data, let $g_{t,s}$,
$t\in\lbrack0,\varepsilon)$, solve the Ricci-Hamilton-DeTurck flow
$\frac{\partial}{\partial t}g_{t,s}=-2\operatorname{Ric}{}_{g_{t,s}
}+\mathcal{L}_{W_{t,s}}g_{t,s}$, where $W_{t,s}=\operatorname{tr}^{1,2}
{}_{g_{t,s}}(\nabla_{g_{t,s}}-\nabla_{g_{t,0}})$. Note that $\frac{\partial
}{\partial t}g_{t,0}=-2\operatorname{Ric}{}_{g_{t,0}}$. Let $v_{t,s}
=\frac{\partial}{\partial s}g_{t,s}$. We have $\frac{\partial^{2}}{\partial
s\partial t}|_{s=0}g_{t,s}=\Delta_{L}v_{t,0}$, which equals $\frac
{\partial^{2}}{\partial t\partial s}|_{s=0}g_{t,s}=\frac{\partial}{\partial
t}v_{t,0}$ (Lichnerowicz Laplacian heat equation). We compute $\frac{\partial
}{\partial t}\operatorname{Vol}(g_{t,s})=\frac{1}{2}\int\operatorname{tr}
{}_{g_{t,s}}(\frac{\partial}{\partial t}g_{t,s})d\mu_{g_{t,s}}=-\int
R_{g_{t,s}}d\mu_{g_{t,s}}$ since $\int\operatorname{tr}{}_{g_{t,s}
}(\mathcal{L}_{W_{t,s}}g_{t,s})d\mu_{g_{t,s}}=0$ by the divergence theorem.
Now
\begin{align*}
\frac{\partial}{\partial s}|_{s=0}\int R_{g_{t,s}}d\mu_{g_{t,s}}  &
=-\frac{\partial^{2}}{\partial s\partial t}|_{s=0}\operatorname{Vol}
(g_{t,s})=-\frac{\partial^{2}}{\partial t\partial s}|_{s=0}\operatorname{Vol}
(g_{t,s})\\
& =-\frac{1}{2}\frac{\partial}{\partial t}\int\operatorname{tr}{}_{g_{t,0}
}(v_{t,0})d\mu_{g_{t,0}}\\
& =\int\langle-\operatorname{Ric}{}_{g_{t,0}}+\frac{g_{t,0}}{2}R_{g_{t,0}
},v_{t,0}\rangle_{g_{t,0}}d\mu_{g_{g_{t,0}}}
\end{align*}
since $\int\operatorname{tr}{}_{g_{t,0}}(\frac{\partial}{\partial t}
v_{t,0})d\mu_{g_{t,0}}=\int\operatorname{tr}{}_{g_{t,0}}(\Delta_{L}
v_{t,0})d\mu_{g_{t,0}}=\int\Delta_{g_{t,0}}(\operatorname{tr}{}_{g_{t,0}
}(v_{t,0}))d\mu_{g_{t,0}}=0$. Finally, take $t=0$.
December 18, 2013. The notion of volume, in various guises, occurs throughout
the study of Ricci flow, especially in Perelman's work. Now, per unit increase
in scale $t$, the volume form of a metric changes with velocity $\frac
{\partial}{\partial t}d\mu=-Rd\mu$. By Clairaut's theorem, the variation of
$-Rd\mu$ is equal to the change per unit increase in scale of the variation of
the volume form, i.e.,
$$
\frac{\partial}{\partial t}(\frac{\operatorname{tr}_{g}v}{2}d\mu_{g})=\left(
\langle\operatorname{Ric}-\frac{1}{2}Rg,v\rangle+\operatorname{div}
(\frac{\nabla\operatorname{tr}v}{2})\right)  d\mu.
$$
In the $f$-warped or entropy version of this, we have $\frac{\partial
}{\partial t}(fe^{-f}d\mu)=(-R-\Delta f)e^{-f}d\mu$ under $\frac{\partial
}{\partial t}g=-2(\operatorname{Ric}+\nabla^{2}f)$ and $\frac{\partial
f}{\partial t}=-R-\Delta f$. Integrating this yields that $\mathcal{N}
\doteqdot\int_{\mathcal{M}}fe^{-f}d\mu$ satisfies $-\frac{d\mathcal{N}}
{dt}=\mathcal{F}\doteqdot\int(R+|\nabla f|^{2})e^{-f}d\mu$. If $\frac
{\partial}{\partial s}g=v$ and $\frac{\partial f}{\partial s}=\frac
{\operatorname{tr}_{g}v}{2}$, then the variation of the energy integrand is
\begin{align*}
& \frac{\partial}{\partial s}((-R-\Delta f)e^{-f}d\mu)\\
& =\left(  (-L(v,\nabla f)+2\langle\operatorname{Ric}+\nabla^{2}
f,v\rangle)e^{-f}+\operatorname{div}(e^{-f}\{\frac{\nabla\operatorname{tr}
v}{2}-v(\nabla f)\})\right)  d\mu\doteqdot A,
\end{align*}
where $L(v,X)=\operatorname{div}^{2}v+\langle\operatorname{Ric},v\rangle
-2\langle\operatorname{div}v,X\rangle+v(X,X)$ is the linear trace Harnack
quadratic. On the other hand, $\frac{\partial}{\partial s}(fe^{-f}d\mu
)=\frac{\operatorname{tr}_{g}v}{2}e^{-f}d\mu$. So Perelman's version is
$\frac{\partial}{\partial t}(\frac{\operatorname{tr}_{g}v}{2}e^{-f}d\mu
)=\frac{\partial^{2}}{\partial s\partial t}(fe^{-f}d\mu)=A$, using Clairaut's
theorem. Note that integration by parts gives $\int L(v,\nabla f)e^{-f}
d\mu=\int\langle\operatorname{Ric}+\nabla^{2}f,v\rangle e^{-f}d\mu$, from
which one obtains Perelman's energy variation formula.
In Section 6.2 of arXiv:0211159 Perelman argues that the $\mathcal{W}$-entropy
(i.e., $\mathcal{F}$ with scaling) integrand is a warped scalar curvature. So,
without scaling (i.e., $\tau$), we would have the correspondences
$\mathcal{N}\sim\operatorname{Vol}$ and $\mathcal{F}\sim\int Rd\mu$, which is
also clear from taking $f=\operatorname{const}\neq0$ as a special case.
However, in 6.2, Perelman's volume is essentially $\int e^{-f}d\mu$, which is
constant under the above variations.
