I'm taking the variation of the Gauss-Bonnet action $$\mathcal{L}_{GB} = \frac{1}{2}\left(R^{2} - 4R^{\mu\nu}R_{\mu\nu} + R^{\mu\nu\alpha\beta}R_{\mu\nu\alpha\beta}\right)$$ to obtain the equations of motion and the boundary terms. The Equations of motion (E.O.M.) I've end up with are those found in the literature $$\left[2(RR_{\mu\nu} - 2R^{\alpha\beta}R_{\mu\alpha\nu\beta} - 2R_{\mu}{}^{\alpha}R_{\nu\alpha} + R_{\mu\rho\alpha\beta}R_{\nu}{}^{\rho\alpha\beta}) - \frac{1}{2}g_{\mu\nu}\mathcal{L}_{GB}\right]\delta g^{\mu\nu}.$$ But I'm not sure how to deal with the remaining terms $$\left(Rg^{\mu\alpha}g^{\nu\beta} - 4R^{\nu\beta}g^{\mu\alpha} + R^{\mu\nu\alpha\beta}\right)\delta R_{\mu\nu\alpha\beta}.$$ My first tought was to make this a total derivative to make it a boundary term (which seems reasonable since I haven't found any article including them in the E.O.M.)... But I can't find a "Palatini Identity" for $\delta R_{\mu\nu\alpha\beta}$.
By this I mean that since $$ \delta R_{b c d}^{a}=\nabla_{c}\left(\delta \Gamma_{b d}^{a}\right)-\nabla_{d}\left(\delta \Gamma_{b c}^{a}\right) $$ $$\delta R_{b d}=\nabla_{a}\left(\delta \Gamma_{b d}^{a}\right)-\nabla_{d}\left(\delta \Gamma_{b a}^{a}\right)$$ there might be an equation like $$\delta R_{\rho\sigma\mu\nu} \overset{?}{=} \nabla_{\mu}\left(\delta\Gamma_{\rho\nu\sigma}\right) - \nabla_{\nu}\left(\delta\Gamma_{\rho\mu\sigma}\right).$$ But I haven't been able to derive such equation. Is there a type of "Palatini identity" for the purely covariant Riemann tensor or any expression that can help me make these terms a total derivative?
