$$\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$$