Deriving Gauss-Bonnet Gravity (Or just higher order corrections)

You can derive the desired expression in the following way:

\begin{align}\delta(R_{ab}R^{ab}) &=\delta R_{ab} R^{ab}+R_{ab}\delta R^{ab}\\ &=\delta R_{ab}R^{ab}+R_{ab}\delta R_{cd}g^{ca}g^{db}\\ &=\delta R_{ab}R^{ab}+R^{cd}\delta R_{cd}\\ &=2R^{ab}\delta R_{ab}\\ &=2R^{ab}\delta(R_{cadb}g^{cd})\\ &=2R^{ab}\delta(R_{acbd}g^{cd})\qquad\mbox{(by symmetry of the Riemann tensor)}\\ &=2R^{ab}\delta R_{acbd}g^{cd}+2R^{ab}R_{acbd}\delta g^{cd}. \end{align} In the last line, the variation of the Riemann tensor vanishes due to the assumption that it is constant, leaving us with

$$\delta(R_{ab}R^{ab})=2R^{ab}R_{acbd}\delta g^{cd}.$$

The way to resolving this is as follows. First, don't assume that the Riemann tensor is a constant:

\begin{align*} \delta(R_{\mu\nu}R^{\mu\nu}) & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha} + \square R_{\mu\nu} -\nabla_{\alpha}\nabla_{\mu}R_{\nu}{}^\alpha - \nabla_{\alpha}\nabla_{\nu}R_{\mu}{}^\alpha + g_{\mu\nu} \nabla_{\alpha}\nabla_{\beta}R^{\alpha\beta}]\delta g^{\mu\nu} \\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha}-((\nabla_{\mu}\nabla_{\alpha}+ \square R_{\mu\nu} -[\nabla_\mu,\nabla_\alpha])R_{\nu}{}^\alpha \\&\quad+(\nabla_{\nu}\nabla_{\alpha}-[\nabla_\nu,\nabla_\alpha])R_{\mu}{}^\alpha) + g_{\mu\nu} \nabla_{\alpha}\nabla_{\beta}R^{\alpha\beta}]\delta g^{\mu\nu}\\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha}+ \square R_{\mu\nu} \\&\quad-\nabla_\mu\nabla_\nu R + [\nabla_\mu,\nabla_\alpha]R_{\nu}{}^\alpha + [\nabla_\nu,\nabla_\alpha]R_{\mu}{}^\alpha + \frac{1}{2}g_{\mu\nu} \square R]\delta g^{\mu\nu}\\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha}+ \square R_{\mu\nu} -\nabla_\mu\nabla_\nu R +2(-R_{\alpha \nu}R^{\alpha}{}_\mu \\& \quad+ R^{\alpha\beta}R_{\mu\alpha\nu\beta})+ \frac{1}{2}g_{\mu\nu} \square R]\delta g^{\mu\nu} \\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} -\nabla_\mu\nabla_\nu R + \square R_{\mu\nu} +2R^{\alpha\beta}R_{\mu\alpha\nu\beta}+ \frac{1}{2}g_{\mu\nu} \square R]\delta g^{\mu\nu}. \end{align*}

Dropping the constant terms, we arrive at $ -\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta}+2R^{\alpha\beta}R_{\mu\alpha\nu\beta}$, including the term from varying $\sqrt{|g|}$. The problem with the previous solution is in the transition from the second to the third line, because one still has to vary the metric, and one should vary everything first without assuming constants, since the commutator of the covariant derivatives do not commute. I think the problem arises because even if the Riemann tensor is a constant, the metric that comes up with it does not necessarily (e.g. $S^2$) that it is constant, so we can't just say $\delta R_{abcd} = 0$, and we have to check the parts of it.