  We report constraints on the three-dimensional orbital architecture for all
four planets known to orbit the nearby M dwarf Gliese 876 based solely on
Doppler measurements and demanding long-term orbital stability. Our dataset
incorporates publicly available radial velocities taken with the ELODIE and
CORALIE spectrographs, HARPS, and Keck HIRES as well as previously unpublished
HIRES velocities. We first quantitatively assess the validity of the planets
thought to orbit GJ 876 by computing the Bayes factors for a variety of
different coplanar models using an importance sampling algorithm. We find that
a four-planet model is preferred over a three-planet model. Next, we apply a
Newtonian MCMC algorithm to perform a Bayesian analysis of the planet masses
and orbits using an n-body model in three-dimensional space. Based on the
radial velocities alone, we find that a 99% credible interval provides upper
limits on the mutual inclinations for the three resonant planets
($\Phi_{cb}<6.20^\circ$ for the "c" and "b" pair and $\Phi_{be}<28.5^\circ$ for
the "b" and "e" pair). Subsequent dynamical integrations of our posterior
sample find that the GJ 876 planets must be roughly coplanar
($\Phi_{cb}<2.60^\circ$ and $\Phi_{be}<7.87^\circ$), suggesting the amount of
planet-planet scattering in the system has been low. We investigate the
distribution of the respective resonant arguments of each planet pair and find
that at least one argument for each planet pair and the Laplace argument
librate. The libration amplitudes in our three-dimensional orbital model
supports the idea of the outer-three planets having undergone significant past
disk migration.
