  We develop a Monte Carlo Comptonization model for the X-ray spectrum of
accretion-powered pulsars. Simple, spherical, thermal Comptonization models
give harder spectra for higher optical depth, while the observational data from
Vela X-1 show that the spectra are harder at higher luminosity. This suggests a
physical interpretation where the optical depth of the accreting plasma
increases with mass accretion rate. We develop a detailed Monte-Carlo model of
the accretion flow, including the effects of the strong magnetic field ($\sim
10^{12}$ G) both in geometrically constraining the flow into an accretion
column, and in reducing the cross section. We treat bulk-motion Comptonization
of the infalling material as well as thermal Comptonization. These model
spectra can match the observed broad-band {\it Suzaku} data from Vela X-1 over
a wide range of mass accretion rates. The model can also explain the so-called
"low state", in which the uminosity decreases by an order of magnitude. Here,
thermal Comptonization should be negligible, so the spectrum instead is
dominated by bulk-motion Comptonization.
