  We develop and analyse finite volume methods for the Poisson problem with
boundary conditions involving oblique derivatives. We design a generic
framework, for finite volume discretisations of such models, in which internal
fluxes are not assumed to have a specific form, but only to satisfy some
(usual) coercivity and consistency properties. The oblique boundary conditions
are split into a normal component, which directly appears in the flux balance
on control volumes touching the domain boundary, and a tangential component
which is managed as an advection term on the boundary. This advection term is
discretised using a finite volume method based on a centred discretisation (to
ensure optimal rates of convergence) and stabilised using a vanishing boundary
viscosity. A convergence analysis, based on the 3rd Strang Lemma \cite{DPD18},
is conducted in this generic finite volume framework, and yields the expected
$\mathcal O(h)$ optimal convergence rate in discrete energy norm.
  We then describe a specific choice of numerical fluxes, based on a
generalised hexahedral meshing of the computational domain. These fluxes are a
corrected version of fluxes originally introduced in \cite{Medla.et.al2018}. We
identify mesh regularity parameters that ensure that these fluxes satisfy the
required coercivity and consistency properties. The theoretical rates of
convergence are illustrated by an extensive set of 3D numerical tests,
including some conducted with two variants of the proposed scheme. A test
involving real-world data measuring the disturbing potential in Earth gravity
modelling over Slovakia is also presented.
