The {100} and (110) cleavage energies of body-centered cubic TiZrNbHf high-entropy alloy are calculated using two alloy models: special quasi-random structures (SQSs) and the coherent potential approximation (CPA). The projector augmented wave method, as implemented in the Vienna ab initio simulation package (VASP), in combination with SQSs is adopted to evaluate the impact of local lattice distortions, whereas the exact muffin-tin orbitals (EMTO) method is used in combination with both SQSs and CPA to study the effect of chemical disorder using rigid underlying lattices. The variations of the cleavage energy as a function of surface chemistry and structure from the EMTO and VASP calculations are consistent with each other. Furthermore, the cleavage energies from CPA are in good agreement with those from SQSs, confirming that an averaged supercell approach reproduces well the mean-field CPA results. The alloy's cleavage energies estimated by the rule of mixtures compare well with those from the direct calculations, and the surface chemistry dependence of the cleavage energies is mainly controlled by the number of Nb atoms in the surface terminal layers owing to the large cleavage energy of Nb metal.