In this paper, we elaborate on the operator effective medium approximation developed recently in Popov et al. [Phys. Rev. B 94, 085428 ( 2016)] to get insight into the surface polariton excitation at the interface of a multilayer hyperbolic metamaterial (HMM). In particular, we find that HMMs with bilayer unit cells support the TE- and TM-polarized surface waves beyond the Maxwell Garnett approximation due to the spatial dispersion interpreted as effective magnetoelectric coupling. The latter is also responsible for the dependence of surface wave propagation on the order of layers in the unit cell. Elimination of the magnetoelectric coupling in three-layer unit cells complying with inversion symmetry restores the qualitative regularity of the Maxwell Garnett approximation, as well as strongly suppresses the influence of the order of layers in the unit cell.