We apply the Hermitian phase operator method (HPOM), as discussed by Pegg and Barnett [D. T. Pegg and S. M. Barnett, Phys. Rev. A 39, 1665 (1989)] to a two-well Bose-Hubbard model. Application of the HPOM formalism yields an approximate quantum phase model, or a Schr\"odinger-like differential equation for the phase variable $\ensuremath{\theta}$, which is a conjugate variable to the half of the number difference between the wells, $m=({N}_{1}\ensuremath{-}{N}_{2})∕2$. In the construction of the model we take care so that the Hermiticity of the original Bose-Hubbard model in number representation is inherited. We demonstrate that the quantum phase model supersedes two theoretical models suggested for the group of so-called ``normal'' states, for which $⟨\ensuremath{\theta}⟩=⟨m⟩=0$: a non-Hermitian ``exact quantum phase model'' by Anglin, Drummond, and Smerzi [J. R. Anglin, P. Drummond, and A. Smerzi, Phys. Rev. A 64, 063605 (2001)], and the ``phonon model'' by Javanainen and Ivanov [J. Javanainen and M. Yu. Ivanov, Phys. Rev. A 60, 2351 (1999)]. We apply the HPOM to the so called ``superfluid'' states, for which $⟨m⟩=0$ and $⟨\ensuremath{\theta}⟩=\ensuremath{\pi}$, and find their superfluid phonon model. We find that the superfluid phonon model has a threshold point $\mathcal{M}=\ensuremath{\delta}∕(4gN)=1$, at which the model diverges. Here, $\ensuremath{\delta}$ is the tunneling rate between the wells and $g$ is the self-interaction, while $N$ is the total number of particles. The above threshold phonon model is quite accurate in the description of the number and phase statistics of the superfluid states. We show that at the threshold point HPOM itself can be exactly solved. However, the agreement between the exact solution and the HPOM solution is poor: $⟨{m}^{2}⟩$ close to the threshold assumes values $\ensuremath{\sim}{N}^{2}$, thus voiding the second order expansion the HPOM is based on, $⟨{m}^{2}⟩⪡{N}^{2}$.