Suppose we are given the Stokes equations with Neumann conditions on part of the boundary:
$-\nabla\cdot\boldsymbol{\sigma} = \mathbf{f}, \quad \text{and} \quad \nabla\cdot \mathbf{u} = 0 \quad \text{in} \quad \Omega$
Here $\boldsymbol{\sigma} := \nu(\nabla \mathbf{u} + (\nabla \mathbf{u})^T) - p I$, and $\mathbf{u}=[u_1,\dots,u_d]^T$. Additionally, we assume the boundary conditions of Dirichlet and Neumann-type
$\mathbf{u} = 0 \quad \text{on} \quad \partial\Omega\backslash\Gamma, \quad \text{and} \quad \boldsymbol{\sigma}\cdot\mathbf{n} = \mathbf{t}\quad \text{on} \quad \Gamma$,
where $\mathbf{n}$ denotes the outward pointing normal vector. Assuming that $\nu$ is constant in the whole of $\Omega$, we can rewrite the momentum equation as
$-\nu\Delta\mathbf{u} + \nabla p = \mathbf{f}, \quad \text{and} \quad \nabla\cdot \mathbf{u} = 0 \quad \text{in} \quad \Omega$
by using the incompressibility $\nabla \cdot \mathbf{u} = 0$ and some elementary manipulations.
What puzzles me now is the following:
By the usual procedure, we obtain the following weak form: find $(\mathbf{u},p) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$ such that
$\int_\Omega (\nu\nabla\mathbf{u} : \nabla\mathbf{v} - p\nabla\cdot\mathbf{v} - q\nabla\cdot\mathbf{u}) \,{\rm d}x - \int_\Gamma (\nu\nabla\mathbf{u}-pI)\cdot\mathbf{n}\cdot\mathbf{v}\,{\rm d}s = \int_\Omega \mathbf{f}\cdot\mathbf{v}\,{\rm d}x$
for all test functions $(\mathbf{v},q) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$. At this point, one would normally replace the remaining boundary terms by the Neumann boundary conditions and put them to the right hand side. However, when considering the physical traction boundary conditions given above, we clearly have that
$\mathbf{t} = \boldsymbol{\sigma}\cdot\mathbf{n} = (\nu(\nabla \mathbf{u} + (\nabla \mathbf{u})^T)\cdot\mathbf{n} - p\mathbf{n} \ne \nu\nabla \mathbf{u}\cdot\mathbf{n} - p\mathbf{n}$
unless for some reason $(\nabla \mathbf{u})^T\cdot\mathbf{n} = 0$ of course. This now leads to the weak form: find $(\mathbf{u},p) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$ such that
$\int_\Omega (\nu\nabla\mathbf{u} : \nabla\mathbf{v} - p\nabla\cdot\mathbf{v} - q\nabla\cdot\mathbf{u}) \,{\rm d}x + \int_\Gamma \nu(\nabla\mathbf{u})^T\cdot\mathbf{n}\cdot\mathbf{v}\,{\rm d}s = \int_\Omega \mathbf{f}\cdot\mathbf{v}\,{\rm d}x - \int_\Gamma \mathbf{t}\cdot\mathbf{v}\,{\rm d}s$.
for all test functions $(\mathbf{v},q) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$. However, what I am really interested in is the problem in finite dimensional spaces $\mathbf{V}_h \times Q_h \subset \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$ for the velocity and pressure, respectively. Let us assume that these are LBB stable. Following the standard procedure for showing stability, one would still have to show that:
$a(\mathbf{u},\mathbf{v}) := \int_\Omega \nu\nabla\mathbf{u} : \nabla\mathbf{v} \,{\rm d}x + \int_\Gamma \nu (\nabla\mathbf{u})^T\cdot\mathbf{n}\cdot\mathbf{v}\,{\rm d}s$.
is uniformly coercive on the kernel space $\mathbf{\tilde V}_h := \{ \mathbf{v}_h \in \mathbf{V}_h : b(\mathbf{v}_h,q_h) = 0, \forall q\in Q_h \}.$ This is something, which I could not find in the literature, and did not manage to show myself. Perhaps some stabilization terms or sophisticated analysis techniques might be needed.
If one wants to discretize the above problem by the finite element method, one can of course always discretize the alternative form
${\tilde a}(\mathbf{u},\mathbf{v}) := \int_\Omega \tfrac{\nu}{2}(\nabla\mathbf{u} + (\nabla\mathbf{u})^T) : (\nabla\mathbf{v} + (\nabla\mathbf{v})^T)\,{\rm d}x$
where the traction boundaries are no problem to incorporate during the derivation.
However, for constant viscosities it might sometimes be advantageous to not implement the above form with the cross-derivatives involved. Reasons may be that the effort to implement this in some legacy code might be too large or efficiency concerns due to the additional couplings, etc.
tl;dr: Most people who use the Laplacian-Stokes only consider pseudo-traction boundary conditions, i.e., they impose $\nu\nabla \mathbf{u}\cdot\mathbf{n} - p\mathbf{n} = \mathbf{t}$, where it assumed that $(\nabla \mathbf{u})^T\cdot\mathbf{n} \approx 0$. I kind of doubt that there is always a physical justification for this, and rather expect that this is mostly done for convenience. However, I would be surprised if nobody ever investigated the use of the Laplacian form with full traction conditions. I was not very successful yet with my literature study and hope that somebody might have additional resources or some more information for me.
What I have found so far is:
Laplace form of Navier-Stokes equations: A safe path or wrong way? by A. Limache and S. Idelsohn (2006)
This paper raises some similar questions.