This is really just a long commentary about your question. First, it is always possible to write everything without using coordinates, because the indices can refer to a (moving) frame of tangent vectors. If I understand correctly, your main goal is to not use indices.
I've always preferred index-free formulas over ones with indices. But usually only for the final formula. It is unusual in Riemannian geometry to be able to write full calculations without indices.
If you want to compute the variation of something, call it $\mathrm{Blah}$, with respect to something else, say $\mathrm{Else}$, then you need a precise definition of $\mathrm{Blah}$ with respect to $\mathrm{Else}$. Gauge theory is in a sense less nonlinear than Riemannian geometry. So it is often possible to write rigorous formulas without indices and that can be differentiated with respect to a variation in the connection relatively easily.
However, this simply isn't true in Riemannian geometry, and your two examples demonstrate this well. The concepts of determinant and trace (with respect to a Riemannian metric) are awkward to define rigorously using index-free notation, i.e., without using a basis of tangent vectors.
So before you can even differentiate the volume form without indices, you need a formula for it that does not use indices.
Another awkward issue is when you want to contract a tensor of higher order with a lower order tensor. Compare $$ g^{jk}\nabla_jR_{kl}\,dx^l $$ to $$ g^{jk}\nabla_lR_{jk}\,dx^l. $$ Each is a contraction of the tensors $g^{-1}$ and $\nabla R$. So what notation can you use do distinguish between these two possibilities without using indices?
It might be possible to invent notation that allows you to do calculation without using indices, but I don't know of any successful effort to do this. The closest I know of is Penrose's abstract index notation.
$\newcommand\Hom{\operatorname{Hom}}$ ADDED: You can of course define everything functorially.
In particular, if $T = T_xM$, then $g \in S^2T^*$ defines functorially a map $g: T \rightarrow T^*$, which induces functorially a map $$\det g: \Lambda^nT \rightarrow \Lambda^nT^*.$$ This means that $\det g \in \Lambda^nT^*\otimes\Lambda^nT^*$. If $g$ is positive definite (i.e., $g(v,v) > 0$ for any $v \ne 0$), then for any nonzero $\omega \in \Lambda^nT$, $$(\det g)(\omega,\omega) > 0.$$ Since $\dim \Lambda^nT^* = 1$, this implies that there exists, unique up to sign, $\omega \in \Lambda^nT^*$ such that $$\omega\otimes \omega = \det g.$$ Then $$dV_g = \omega.$$ I'm pretty sure you could do all the variation calculations using this definition of $dV_g$, but, as far as I can tell, it isn't worth the trouble.