general relativity - Variation of the purely covariant Riemann tensor

I need to find the variation of the purely covariant Riemann tensor with respect to the metric $g^{\mu \nu}$, i.e. $\delta R_{\rho \sigma \mu \nu}$.

I know that, $R_{\rho \sigma \mu \nu} = g_{\rho \eta} R^\eta_{\sigma \mu \nu}$

and $\delta R^\eta_{\sigma \mu \nu} = \nabla_\mu (\delta \Gamma^\eta_{\nu \sigma})- \nabla_\nu (\delta \Gamma^{\eta}_{\mu \sigma}).$

So do I need to use the product rule like so, $\delta R_{\rho \sigma \mu \nu} = \delta (g_{\rho \eta} R^\eta_{\sigma \mu \nu}) = \delta (g_{\rho \eta}) R^\eta_{\sigma \mu \nu} + g_{\rho \eta} \delta( R^\eta_{\sigma \mu \nu})$,

which gives me a relatively messy answer. Is there another way to approach this?