Modeling results of volcanic plumes on Jupiter's moon Io are presented. Two types of low density axisymmetric SO 2 plume flows are modeled using the direct simulation Monte Carlo (DSMC) method. Thermal radiation from all three vibrational bands and overall rotational lines of SO 2 molecules is modeled. A high resolution computation of the flow in the vicinity of the vent was obtained by multidomain sequential calculation to improve the modeling of the radiation signature. The radiation features are examined both by calculating infrared emission spectra along different lines-of-sight through the plume and with the DSMC modeled emission images of the whole flow field. It is found that most of the radiation originates in the vicinity of the vent, and non-LTE (non-local-thermodynamic equilibrium) cooling by SO 2 rotation lines exceeds cooling in the v 2 vibrational band at high altitude.In addition to the general shape of the plumes, the calculated average SO 2 column density (~10 1 6 cm - 2 ) over a Pele-type plume and the related frost-deposition ring structure (at R ~ 500 km from the vent) are in agreement with observations. These comparisons partially validate the modeling. It is suggested that an observation with spatial resolution of less than 30 km is needed to measure the large spatial variation of SO 2 near a Pele-type plume center. It is also found that an influx of 1.1 x 10 2 9 SO 2 s - 1 (or 1.1 x 10 4 kg s - 1 ) is sufficient to reproduce the observed SO 2 column density at Pele. The simulation results also show some interesting features such as a multiple bounce shock structure around Prometheus-type plumes and the frost depletion by plume-induced erosion on the sunlit side of Io. The model predicts the existence of a canopy shock, a ballistic region inside the Pele-type plume, and the negligible effect of surface heating by plume emission.