  In the collisional region at finite temperatures, the collective modes of
superfluids are described by the Landau two-fluid hydrodynamic equations. This
region can now be probed over the entire BCS-BEC crossover in trapped Fermi
superfluids with a Feshbach resonance, including the unitarity region. Building
on the approach initiated by Zaremba, Nikuni, and Griffin in 1999 for trapped
atomic Bose gases, we present a new variational formulation of two-fluid
hydrodynamic collective modes based on the work of Zilsel in 1950 developed for
superfluid helium. Assuming a simple variational ansatz for the superfluid and
normal fluid velocities, the frequencies of the hydrodynamic modes are given by
solutions of coupled algebraic equations, with constants only involving spatial
integrals over various equilibrium thermodynamic derivatives. This variational
approach is both simpler and more physical than a direct attempt to solve the
Landau two-fluid differential equations. Our two-fluid results are shown to
reduce to those of Pitaevskii and Stringari for a pure superfluid at T=0.
