$\operatorname{div} (X \otimes Y)$ using only geometric calculus. It is easy to compute $\operatorname{div} (X \wedge Y)$ using definition $\operatorname{div} X = (\delta (X^\flat))^\sharp$ and primitive identities which consists signs $\star, d, \delta, \flat, \wedge$. But can we compute $\operatorname{div} (X \otimes Y)$ in such way? It's seems much more difficult because we can't define codifferential $\delta$ on any $k$-form from $\Gamma^\infty(T^{*k} M)$.
 A: This isn't a complete answer, I'm not sure how to finish off/simplify the argument below
Recall that we want the divergence of a tensor field to simply be the negative adjoint of the covariant derivative, $\nabla$, with respect to the integrated inner product. This simply means that $\operatorname{div}=-\nabla^*$ where $\nabla^*$ satisfies the following: if $S$ is a compactly supported $(r,s)$-type tensor field and $T$ is a compactly supported $(r,s+1)$-type tensor field then
$$\int_{M}\langle\nabla S,T\rangle~\operatorname{vol}_g=\int_{M}\langle S,\nabla^*T\rangle~\operatorname{vol}_g.$$
It isn't hard to check that if $\{E_1,\dots,E_n\}$ is an orthonormal frame then
$$(\operatorname{div}(T))(Z_1,\dots,Z_{s-1},\alpha_1,\dots,\alpha_r)=\sum(\nabla_{E_i}T)(E_i,Z_1,\dots,Z_{s-1},\alpha_1,\dots,\alpha_r).$$
Now we consider the case when $T=X\otimes Y$, $X\in T^{s_1}_{r_1}M$ and $Y\in T^{s_2}_{r_2}M$. For the sake of conciseness, for any tuple $\vec{v}=(v_1,\dots,v_\ell)$, we write $\vec{v}_{i,j}=(v_i,\dots,v_j)$. Let $\vec{Z}=(Z_1,\dots,Z_{s-1})$ be a tuple of vector fields, and $\vec{\alpha}=(\alpha_1,\dots,\alpha_r)$ be a tuple of 1-forms. 
We compute 
\begin{align}
(\operatorname{div}(X\otimes Y))(\vec{Z}_{1,s-1},\vec{\alpha}_{1,r}) & =\sum(\nabla_{E_i}X\otimes Y + X\otimes\nabla_{E_i}Y)(E_i,\vec{Z}_{1,s-1},\vec{\alpha}_{1,r}) \\
& = \sum (\nabla_{E_i}X)(E_i,\vec{Z}_{1,s_1-1},\vec{\alpha}_{1,r_1})Y(\vec{Z}_{s_1,s-1},\vec{\alpha}_{r_1+1,r})\\
& \qquad +X(E_i,\vec{Z}_{1,s_1-1},\vec{\alpha}_{1,r})(\nabla_{E_i}Y)(\vec{Z}_{s_1,s},\vec{\alpha}_{r_1+1,r})\\
& = \operatorname{div}(X)(\vec{Z}_{1,s_1-1},\vec{\alpha}_{1,r_1})Y(Z_{s_1,s-1},\vec{\alpha}_{r_1+1,r})\\
& \qquad + \sum X(E_i,\vec{Z}_{1,s_1-1},\vec{\alpha}_{1,r_1})(\nabla_{E_i}Y)(\vec{Z}_{s_1,s},\vec{\alpha}_{r_1+1,r}) \\
& = (\operatorname{div}(X)\otimes Y)(Z_{1,s-1},\vec{\alpha}_{1,r})\\
& \qquad \color{red}{+ \sum X(E_i,\vec{Z}_{1,s_1-1},\vec{\alpha}_{1,r_1})(\nabla_{E_i}Y)(\vec{Z}_{s_1,s},\vec{\alpha}_{r_1+1,r})}
\end{align}
I'm not sure how to simplify the term above in red in such a way that we don't need to use the natural isomorphism of tensor fields with $\mathscr{C}^\infty(M)$-multilinear maps. I hope this helps though! Maybe somebody else will see how to simplify the red term I have above.
