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.