What is wrong with the derivation? Let $(M^n,g)$ be a Riemannian manifold, and $T$ a symmetric $(1,1)$-tensor field, i.e., $\langle T(X),Y\rangle = \langle X,T(Y)\rangle $. For convenience, denote $$\Delta_Tu=\sum_i\langle \nabla_{e_i}\nabla u, Te_i\rangle $$
and
$$\mathrm{Ric}_T(X,Y)=\sum_i\langle R(X,e_i)(Te_i), Y\rangle , $$
where $u$ is a smooth function on $M$ and $\{e_i\}$ is a local ON frame field.
Now assume $T$ is a Codazzi operator, i.e., for any $X,Y\in \Gamma(TM)$, $(\nabla_XT)Y=(\nabla_YT)X$. We choose $\{e_i\}_{i=1}^n$ be a local orthonormal frame field of $M$ such that $\nabla_{\star }e_i=0$ at the considered point. For the distance function r(x) from a fixed point $x_0$, by the definition, we have ($\nabla_XT$ is symmetric since $T$ is symmetric)
  \begin{equation*}
    \begin{split}
      \Delta_{\nabla_{\partial_r}T}r=&\sum_{i=1}^n\langle \nabla_{e_i}\partial_r,(\nabla_{\partial_r}T)e_i\rangle=\sum_{i=1}^n\langle \nabla_{e_i}\partial_r,(\nabla_{e_i}T)\partial_r\rangle \\
      =&\sum_{i=1}^ne_i\langle \partial_r,(\nabla_{e_i}T)\partial_r\rangle -\sum_{i=1}^n\langle \partial_r,(\nabla_{e_i}\nabla_{e_i}T)\partial_r\rangle -\sum_{i=1}^n\langle \partial_r,(\nabla_{e_i}T)(\nabla_{e_i}\partial_r)\rangle .
    \end{split}
  \end{equation*}
  However, 
  \begin{equation*}
    \begin{split}
      \sum_{i=1}^n\langle \partial_r,(\nabla_{e_i}T)(\nabla_{e_i}\partial_r)\rangle =&\sum_{i=1}^n\langle (\nabla_{e_i}T)\partial_r,\nabla_{e_i}\partial_r\rangle \\
      =&\sum_{i=1}^n\langle (\nabla_{\partial_r}T)e_i,\nabla_{e_i}\partial_r\rangle =\Delta_{\nabla_{\partial_r}T}r.
    \end{split}
  \end{equation*}
  Hence, we obtain 
  \begin{equation}
    \begin{split}
      \Delta_{\nabla_{\partial_r}T}r=\frac{1}{2}\sum_{i=1}^ne_i\langle \partial_r,(\nabla_{e_i}T)\partial_r\rangle -\frac{1}{2}\sum_{i=1}^n\langle \partial_r,(\nabla_{e_i}\nabla_{e_i}T)\partial_r\rangle 
    \end{split}
  \end{equation}
  We now compute the two terms of the R.H.S. of the above equality. Firstly, notice that $\nabla_{\partial_r}\partial_r=0$, we have
  \begin{equation}
    \begin{split}
      \sum_{i=1}^ne_i\langle \partial_r,(\nabla_{e_i}T)\partial_r\rangle =&\sum_{i=1}^ne_i\langle \partial_r,(\nabla_{\partial_r}T)e_i\rangle =\sum_{i=1}^ne_i\langle (\nabla_{\partial_r}T)\partial_r,e_i\rangle \\
      =&\sum_{i=1}^ne_i (\partial_r\langle T\partial_r, e_i\rangle )-\sum_{i=1}^ne_i\langle T\partial_r,\nabla_{\partial_r}e_i\rangle \\
      =&\sum_{i=1}^n\partial_r(e_i\langle T\partial_r,e_i\rangle )-\sum_{i=1}^n\langle T\partial_r, \nabla_{e_i}\nabla_{\partial_r}e_i\rangle \\ 
      =&\sum_{i=1}^n\partial_r\langle (\nabla_{e_i}T)\partial_r,e_i\rangle +\sum_{i=1}^n\partial_r\langle T\nabla_{e_i}\partial_r, e_i\rangle \\
       &+\sum_{i=1}^n\partial_r\langle T\partial_r,\nabla_{e_i}e_i\rangle -\sum_{i=1}^n\langle T\partial_r, \nabla_{e_i}\nabla_{\partial_r}e_i\rangle \\
      =&\sum_{i=1}^n\langle (\nabla_{\partial_r}\nabla_{e_i}T)\partial_r,e_i\rangle +\partial_r(\Delta_Tr)\\
       &+\sum_{i=1}^n\langle T\partial_r,\nabla_{\partial_r}\nabla_{e_i}e_i\rangle -\sum_{i=1}^n\langle T\partial_r, \nabla_{e_i}\nabla_{\partial_r}e_i\rangle \\
      =&\sum_{i=1}^n\langle (\nabla_{\partial_r}\nabla_{e_i}T)\partial_r,e_i\rangle +\partial_r(\Delta_Tr)+\mathrm{Ric}(\partial_r, T\partial_r).
    \end{split}
  \end{equation}
  Secondly, 
  \begin{equation}
    \begin{split}
      \sum_{i=1}^n\langle \partial_r,(\nabla_{e_i}\nabla_{e_i}T)\partial_r\rangle =&\sum_{i=1}^n\langle \partial_r,\nabla_{e_i}((\nabla_{e_i}T)\partial_r)\rangle -\sum_{i=1}^n\langle \partial_r,(\nabla_{e_i}T)\nabla_{e_i}\partial_r\rangle \\
      =& \sum_{i=1}^n\langle \partial_r,\nabla_{e_i}((\nabla_{\partial_r}T)e_i)\rangle -\sum_{i=1}^n\langle \partial_r,(\nabla_{\nabla_{e_i}\partial_r}T)e_i\rangle \\
      =&\sum_{i=1}^n\langle (\nabla_{e_i}\nabla_{\partial_r}T)e_i-(\nabla_{\nabla_{e_i}\partial_r}T)e_i,\partial_r\rangle \\
      =&\sum_{i=1}^n\langle (\nabla_{\partial_r}\nabla_{e_i}T)e_i,\partial_r\rangle -\sum_{i=1}^n\langle (R(\partial_r,e_i)T)e_i,\partial_r\rangle \\
      =&\sum_{i=1}^n\langle (\nabla_{\partial_r}\nabla_{e_i}T)e_i,\partial_r\rangle +\mathrm{Ric}(\partial_r,T\partial_r)-\mathrm{Ric}_T(\partial_r,\partial_r).
    \end{split}
  \end{equation}
  From the above three equalities we obtain
  \begin{equation*}
    \begin{split}
      \Delta_{\nabla_{\partial_r}T}r
      =\frac{1}{2}\partial_r(\Delta_Tr)
       +\frac{1}{2}\mathrm{Ric}_T(\partial_r,\partial_r).
    \end{split}
  \end{equation*} 
Now, my question is that when $T=\mathrm{Id}_{TM}$ the above equation becomes
 \begin{equation*}
    \begin{split}
    \partial_r(\Delta_r)+\mathrm{Ric}(\partial_r,\partial_r)=0.
    \end{split}
  \end{equation*} 
But it is well known that the Bochner formula for the distance function 
\begin{equation*}
    \begin{split}
    |\mathrm{Hess}r|^2+\partial_r(\Delta_r)+\mathrm{Ric}(\partial_r,\partial_r)=0.
    \end{split}
  \end{equation*} 
This obtain a contradiction. 
What is wrong with the above derivation? Thanks in advence.
 A: Let me re-do the computation using index notation. 
Your goal is to compute $\partial_r (\Delta_T r)$ which we can rewrite as
$$ \nabla^c( T^{ab} \nabla^2_{ab} r) \nabla_c r = \nabla^c T^{ab} \nabla^2_{ab} r \nabla_c r +  T^{ab} (\nabla^c \nabla_a  \nabla_b r) \nabla_c r$$
The first term on the right is $\Delta_{\nabla_{\partial_r} T} r$. So we focus our attention on the second term. 
$$ \nabla^c\nabla_a \nabla_b r = \nabla_a\nabla^c \nabla_b r + [\nabla^c, \nabla_a ] \nabla_b r $$
So (the $\pm$ is just me forgetting which sign is the right one)
$$ T^{ab} (\nabla^c\nabla_a \nabla_b r) \nabla_c r = T^{ab}(\nabla_a \nabla^c \nabla_b r) \nabla_c r \pm \mathrm{Ric}_T(\partial_r, \partial_r) $$
The term in the middle is
$$ \nabla_a \nabla^c \nabla_b r \nabla_c r = \nabla_a (\underbrace{\nabla^c \nabla_b r \cdot \nabla_c r}_{ = 0}) - \nabla^c \nabla_b r \nabla_a \nabla_c r.$$
The term that $=0$ is so due to the fact that $\nabla r$ is geodesic.
So your identity should be
$$ \partial_r (\Delta_T r) + \mathrm{Ric}_T(\partial_r, \partial_r) = \Delta_{\nabla_{\partial_r} T} r - T^{ab} g^{cd} \nabla^2_{bd} r \nabla^2_{ac} r $$
which now reduces correctly to Bochner's identity when $T = g$. 
Remark Notice that nowhere in this derivation is either of the property you started with, namely that $T$ is symmetric and $T$ is Codazzi, used. 
