  We develop a high-order hybridized discontinuous Galerkin (HDG) method for a
linear degenerate elliptic equation arising from a two-phase mixture of mantle
convection or glacier dynamics. We show that the proposed HDG method is
well-posed by using an energy approach. We derive ${\it a priori}$ error
estimates for the proposed HDG method on simplicial meshes in both two- and
three-dimensions. The error analysis shows that the convergence rates are
optimal for both the scaled pressure and the scaled velocity for non-degenerate
problems and are sub-optimal by half order for degenerate ones. Several
numerical results are presented to confirm the theoretical estimates. We also
enhance the HDG solutions by post-processing. The superconvergence rates of
$(k+2)$ and $(k+\frac{3}{2})$ are observed for both a non-degenerate case and a
degenerate case away from the degeneracy. Degenerate problems with low
regularity solutions are also studied, and numerical results show that
high-order methods are beneficial in terms of accuracy.
