  We design, analyse and implement an arbitrary order scheme applicable to
generic meshes for a coupled elliptic-parabolic PDE system describing miscible
displacement in porous media. The discretisation is based on several
adaptations of the Hybrid-High-Order (HHO) method due to Di Pietro et al.
[Computational Methods in Applied Mathematics, 14(4), (2014)]. The equation
governing the pressure is discretised using an adaptation of the HHO method for
variable diffusion, while the discrete concentration equation is based on the
HHO method for advection-diffusion-reaction problems combined with numerically
stable flux reconstructions for the advective velocity that we have derived
using the results of Cockburn et al. [ESAIM: Mathematical Modelling and
Numerical Analysis, 50(3), (2016)]. We perform some rigorous analysis of the
method to demonstrate its $L^2$ stability under the irregular data often
presented by reservoir engineering problems and present several numerical tests
to demonstrate the quality of the results that are produced by the proposed
scheme.
