Zero-gradient boundary condition in a time-independent rotating ocean PDE model

Dear all,

I am aiming to create a code which solves a simple, stationary PDE for sea level. Adopting MATLAB's PDE Toolbox notations, the PDE and coastal boundary conditions for sea level are

which encode respectively vorticity conservation (the physics of the problem) and the constraint of no flow through the coast, with the vector normal to the coastal boundary. The problem is two dimensional (shallow water equations with Coriolis force, see https://en.wikipedia.org/wiki/Shallow_water_equations) and is a 2x2 matrix, . Note and vary horizontally and there is a designated way to take that into account within MATLAB, not shown here for simplicity. The PDE is the 2D "arrested topographic wave" model used in physical oceanography:

with the sea level . The Jacobian term (1st order) dominates the problem and the diffusion term (2nd order) is of much lesser importance.

In my application, it is necessary to have an open boundary to the south (say, ). Such boundaries are "artificial" in oceanography and need to be carefuly thought. Indeed, there is no additional equation constraining the behaviour of the sea level at such an open limit. Nonetheless there is no way to do without the artificial boundary, for many reasons I don't detail here.

The easiest treatment I can think of is to enforce a zero-gradient condition across the southern boundary of my domain, after having arranged that the physics and geometry at proximity of such boundary are simplified. By simplified I mean, any RHS forcing and friction should be nil (respectively and nearby the boundary) and contours perpendicular to the boundary. That way the sea level simply follows contours of and make its way to the boundary (this is trivial from the 1st order method of characteristics) and I could in principle enforce a zero-gradient boundary condition.

The difficulties come from trying to set up a boundary condition in Matlab in this context. The method to control the boundary condition in Matlab does not spit out the gradient I want. That is, it gives me control over (which I don't care about!) instead of . This is because, at the boundary,

The difficulty arise from working with a rotating Earth; and the subsequent orthogonal relationship between pressure-gradient forces and currents in geophysical fluids. Unfortunately we can't stop the Earth from spinning.

So, my question is, how can I enforce a at the boundary when the nearby represented physics () does not involve along the normal direction but rather along the tangential direction?

I hope this question has an answer. I wonder whether adopting a different in the evaluation of the PDE in the domain vs in the evaluation of the gradients at the boundary is a suited method. Alternatively, I wonder whether it is possible to input a tangential boundary condition rather than a normal (Neumann) boundary condition.

Thank you for your help.

Best,

Sam D.

[Torsten]

Not an answer to your problem, but a remark: the shallow water equations are hyperbolic in nature, but the PDE Toolbox is usually applied to parabolic-elliptic systems of equations. Better use a solver suited to solve hyperbolic systems - I recommend CLAWPACK:

https://www.clawpack.org/