Nonlinear Stokes equations for glacier modeling

Question

What is the difference between

$$\operatorname{div}(-\mu(u)\nabla u+pI)=f\tag{1}$$

and

$$\operatorname{div}(-\mu(u)(\nabla u\color{red}{+(\nabla u)^\mathsf{T}})+pI)=f\tag{2}$$

from the physical point of view, in the glacier modeling context?

Does it make sense to model a glacier with the first one?

To explain the notation, $u$ and $p$ denote the velocity and pressure fields, respectively; $\mu(u)$ denotes a velocity-dependent viscosity coefficient (for example, Glen's power flow law). The gradient of $u$ is

$$\nabla u=\left(\begin{array}{cc}\dfrac{\partial u_1}{\partial x}&\dfrac{\partial u_1}{\partial y}\\ \dfrac{\partial u_2}{\partial x}&\dfrac{\partial u_2}{\partial y}\end{array}\right)$$

and $(\nabla u)^\mathsf{T}$ denotes the transpose matrix of $\nabla u$. (If $u$ has 3 variables its gradient is analogous to the 2-variable case). Frequently, $f$ is $f=(0,0,-\rho g)$, where $g$ is the gravitational acceleration constant and $\rho$ the ice density. $I$ denotes the identity matrix and $div$ the divergence of a matrix is the divergence by rows, for example,

$$\operatorname{div}\left(\begin{array}{cc}w_1&w_2\\ s_1& s_2\end{array}\right)=\left(\begin{array}{cc}\dfrac{\partial w_1}{\partial x}+\dfrac{\partial w_2}{\partial y}\\ \dfrac{\partial s_1}{\partial x}+\dfrac{\partial s_2}{\partial y}\end{array}\right).$$

This way, we may rewrote (1) and (2), respectively, as follows:

$$-\operatorname{div}(\mu(u)\nabla u)+\nabla p=f$$

and

$$-\operatorname{div}(\mu(u)(\nabla u\color{red}{+(\nabla u)^\mathsf{T}})+\nabla p=f$$

Answer 1

If your viscosity $\mu(u)$ is constant, then there is no difference between the two formulations because $$ \text{div} (\mu \nabla u^t) = \mu \nabla (\text{div} u) $$ and for incompressible materials $\text{div} u=0$. So the two formulations are identical for constant viscosities. For variable viscosities, this is a bit more complicated and may in fact not be correct.

In general, the second form you show (with the addition of the transpose of the gradient) is the physically correct form because it involves the strain rate $\varepsilon = \frac 12(\nabla u + \nabla u^t)$ which is the physically correct quantity where the gradient $\nabla u$ has no physical reality. (This is because the invariants of $\varepsilon$ are independent of the choice of coordinate system, but the invariants of $\nabla u$ depend on the choice of coordinate system.)

The first form you show is a simplification one can work with if you know that (i) you have an incompressible fluid, and (ii) you have no Neumann type boundary conditions (as the "natural" boundary conditions for the two formulations are different).