  We investigate the compressional modes of cold neutron stars with cores
consisting of superfluid neutrons, superconducting protons and normal fluid
electrons and muons, and crusts that contain superfluid neutrons plus a normal
fluid of (spherical) nuclei and electrons. We develop a two-fluid formalism for
the core that accounts for leptonic buoyancy, and an analogous treatment for
the crust. We adopt the Cowling approximation, neglecting gravitational
perturbations, but include all effects of the background space-time. We
introduce a phenomenological, easily-modified nuclear equation of state which
contains all of the thermodynamic information required to compute the coupled
fluid oscillations, with parameters that are constrained by nuclear physics and
the requirement that the maximum mass of a neutron star is $\geq 2M_{\odot}$.
Using four parametrizations of this equation of state with nuclear
compressibilities $K=230$-$280$ MeV, we calculate the Brunt-V\"{a}is\"{a}l\"{a}
frequency due to leptonic buoyancy, and find the corresponding $g$-mode
frequencies and eigenfunctions. We find that the WKB approximation reproduces
$g$-mode frequencies closely. We examine the dependence of $g$-mode frequencies
on stellar mass, nuclear compressibility and strength of neutron-proton
entrainment, and compare to previous calculations of $g$-mode frequencies due
to leptonic buoyancy. We also compute the $p$-mode spectra, confirming previous
findings that the two fluids behave as if uncoupled except the case of large
entrainment, and show the existence of nearly resonant mode pairs which could
lead to nonlinear $p$-$g$ instabilities even at zero temperature.
