Well, two and a half years later, I stumbled across an answer. The following is a paraphrasing of Petersen's "Riemannian Geometry," page 40.
We work on a Riemannian manifold $(M,\langle \cdot, \cdot \rangle)$ with Levi-Civita connection $\nabla$. Let $\text{Ric}$ denote the symmetric $(1,1)$-tensor $\text{Ric}(v) = \sum R(v,e_i)e_i$, where $\{e_i\}$ is an orthonormal basis. Let $S$ denote the scalar curvature, i.e. $S = \text{tr}(\text{Ric}) = \sum \langle \text{Ric}(e_i), e_i \rangle$. We note that $\nabla S = dS$.
Prop: $dS = 2\,\text{div}(\text{Ric})$.
Proof: Let $\{E_i\}$ be an orthonormal frame at $p \in M$ with $\nabla E_i|_p = 0$. Let $X$ be a vector field with $\nabla X |_p = 0$. From the second Bianchi identity:
$$\begin{align*}
dS(X)(p) = XS(p) & = X \sum \langle \text{Ric}(E_i), E_i \rangle \\
& = X \sum \langle R(E_i, E_j)E_j, E_i\rangle \\
& = \sum \langle \nabla_X [R(E_i, E_j)E_j], E_i \rangle \\
& = \sum \langle (\nabla_X R)(E_i, E_j)E_j, E_i \rangle \\
& = -\sum \langle (\nabla_{E_j}R)(X, E_i)E_j, E_i \rangle - \sum \langle (\nabla_{E_i}R)(E_j,X)E_j, E_i \rangle \\
& = -\sum (\nabla_{E_j}R)(X, E_i, E_j, E_i) - \sum (\nabla_{E_i}R)(E_j, X, E_j, E_i) \\
& = \sum (\nabla_{E_j}R)(E_j, E_i, E_i, X) + \sum (\nabla_{E_j} R)(E_i, E_j, E_j, X) \\
& = 2 \sum (\nabla_{E_j}R)(E_j, E_i, E_i, X) \\
& = 2 \sum \nabla_{E_j}(R(E_j, E_i, E_i, X) \\
& = 2 \sum \nabla_{E_j} \langle \text{Ric}(E_j), X \rangle \\
& = 2 \sum \nabla_{E_j} \langle \text{Ric}(X), E_j \rangle \\
& = 2 \sum \langle \nabla_{E_j}(\text{Ric}(X)), E_j) \rangle \\
& = 2 \sum \langle (\nabla_{E_j}\text{Ric})(X), E_j \rangle \\
& = 2\,\text{div}(\text{Ric})(X)_p.
\end{align*}$$
Best Answer
Computing $|X|^2$ directly one has \begin{align*} |X|^2&=\nabla_iR_{jk}\nabla^iR^{jk}+\frac{1}{9}\nabla_iRg_{jk}\nabla^iRg^{jk}-\frac{2}{3}\nabla_iR_{jk}\nabla^iRg^{jk}\\ &=|\nabla\operatorname{Ric}|^2+\frac{1}{3}|\nabla R|^2-\frac{2}{3}|\nabla R|^2\\ &=|\nabla\operatorname{Ric}|^2-\frac{1}{3}|\nabla R|^2 \end{align*}