  We introduce a new discretization of a mixed formulation of the
incompressible Stokes equations that includes symmetric viscous stresses. The
method is built upon a mass conserving mixed formulation that we recently
studied. The improvement in this work is a new method that directly
approximates the viscous fluid stress $\sigma$, enforcing its symmetry weakly.
The finite element space in which the stress is approximated consists of
matrix-valued functions having continuous "normal-tangential" components across
element interfaces. Stability is achieved by adding certain matrix bubbles that
were introduced earlier in the literature on finite elements for linear
elasticity. Like the earlier work, the new method here approximates the fluid
velocity $u$ using $H(\operatorname{div})$-conforming finite elements, thus
providing exact mass conservation. Our error analysis shows optimal convergence
rates for the pressure and the stress variables. An additional post processing
yields an optimally convergent velocity satisfying exact mass conservation. The
method is also pressure robust.
