Hello Alexis,
The free surface gradient compatibility term is explained in the release notes of version 5.7 (you can download it from this Website). It consists of considering that the advection field has a linear part and a piecewise constant part which comes from the gradient of the free surface. We keep this complicated form of velocity for some time in the algorithm, before resorting to the classical linear discretisation. This helps suppressing inf-sup oscillations (it creates a different discretisation between depth and velocity).
At open ocean boundary conditions, the problem is not well posed if you prescribe only the free surface, and indeed you have sometimes spurious flows coming into the domain (and you cannot complain because you said that the velocity is free, so you get a solution among an infinity). In this case it is better to have an open boundary which is in deep waters (if possible beyond the continental shelf), where the velocity is small. In 3D in one study the problem was solved by adding a small head loss at the boundary. The equivalent in 2D would be adding a local friction (this creates a damping term for spurious velocities).
Another solution is using the Thompson boundary conditions, however sometimes it yields a free surface elevation that is slightly different from the one measured, so people generally stick to the ill-posed problem with prescribed elevations.
With best regards,
Jean-Michel Hervouet