When dealing with the Navier-Stokes equations, one typically models fluids in which the viscous stress depends linearly on the symmetric velocity gradient, $$\tau^{ij}=2\mu(\nabla u)^{(ij)}=\mu\left(\nabla ^iu^j+\nabla^ju^i\right)$$ So when taking the divergence of the stress, as you do when writing the general continuum equations, we get $$\nabla_j\tau^{ij}=\mu\nabla_j\nabla^iu^j+\mu\nabla_j\nabla^j u^i$$ The second term can be recognized as just the vector Laplacian of $\boldsymbol u$. However, the first term, $$\nabla_j\nabla^i u^j$$ Is not so simple. The covariant derivatives $\nabla_j,\nabla^i$ do not necessarily commute. However, if we are working in rectangular coordinates, they do, and one can write $$\nabla_j\nabla^i u^j=\nabla^i\nabla_ju^j=\nabla^i(\operatorname{div}u)=(\operatorname{grad}\operatorname{div}u)^i$$ And hence $$\nabla\cdot \boldsymbol\tau=\mu\big(\nabla(\nabla\cdot\boldsymbol u)+\nabla^2\boldsymbol u\big)$$
So, my question is this: Can we always say $\nabla_j\nabla^i u^j=(\operatorname{grad}\operatorname{div}u)^i$ ? I should hope the answer is yes, because, when dealing with incompressible flows, one uses this identity to discard the $\nabla_j\nabla^i u^j$ term from the divergence of the viscous stress. So, my thinking is, if this term is always zero in Cartesian coordinates, it should also be zero in any other coordinate system. But, this is not clear from the expression.