  Although boundary conditions are mandatory to solve partial differential
equations, they also represent a transfer of information between the domain
being modelled and its surroundings. In the case of isolated or closed systems,
these can be formulated using free- or no-slip conditions. However, for open
systems, the information transferred through the boundaries is essential to the
dynamics of the system and can have a first order impact on its evolution. This
work addresses regional geodynamic modelling simulating the evolution of an
Earth's piece over millions of years by solving non-linear Stokes flow. In this
open system, we introduce a new approach to impose oblique boundary conditions
generalising the Navier-slip boundary conditions to arbitrary directions in 3D.
The method requires defining both slip and stress constraints. The stress
constraint is imposed utilising a coordinate transformation to redefine the
stress tensor along the boundaries according to the arbitrary direction chosen
while for the slip constraint we utilise Nitsche's method in the context of the
finite element method, resulting in a symmetrised and penalised weak form. We
validate our approach through a series of numerical experiments of increasing
complexity, starting with 2D and 3D linear models. Then, we apply those
boundary conditions to a 3D non-linear geodynamic model of oblique extension
that we compare with a standard model utilising Dirichlet boundary conditions.
Our results show that using Dirichlet boundary conditions strongly influences
the evolution of the system and generates artefacts near and along the
boundaries. In comparison, the model using the generalised Navier-slip boundary
conditions behaves closely to a model with an unbounded domain, providing a
physically interpretable solution near and along the boundaries.
