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.
This is a very silly question. Because $\mathbb R^3$ is flat, we can find a coordinate system (e.g, Cartesian) in which the Riemann tensor is zero everywhere. And, due to the following property of the Riemann tensor, $$[\nabla_c,\nabla_d]w_b=w_aR^a{}_{bcd}$$ In the Cartesian coordinate system we get $[\nabla_b,\nabla_c]=0$ as expected. However, when we change from our original coordinate system $\{x^i\}$ to a new coordinate system $\{\bar x^i\}$ the Riemann tensor transforms as $$\overline{R}^a{}_{bcd}=\frac{\partial\bar{x}^a}{\partial x^i}\frac{\partial x^j}{\partial \bar{x}^b}\frac{\partial x^k}{\partial \bar{x}^c}\frac{\partial x^l}{\partial \bar{x}^d}R^i{}_{jkl}$$ The components transform linearly. So if the Riemann tensor is zero in one coordinate system, it has to be zero in any other coordinate system for $\mathbb R^3$ and so we can always commute the covariant derivatives.