This is a great problem.   Getting the boundary conditions right is going to be a real challenge.  They should be able to be implmented as slip on the free boundary, possibly as a weak form contribution. Instead of posting an article, I have attached a document about this sort of problem in 2D from a model implemented in Comsol Multiphysics.  This document provides some basic analysis of the problem, and may be of some use here.