Computational results and comparisons with experimental data are presented for simulations of axisymmetric turbulent argon plasma jets flowing into a cold air environment. The calculations were performed using the LAVA code [J. D. Ramshaw and C. H. Chang, Plasma Chem. Plasma Process. 12, 299 (1992)], and were designed to simulate experiments performed by Brossa and Pfender [Plasma Chem. Plasma Process. 8, 75 (1988)] (BP) and by Fincke et al [private communication, 1992] (FSH). To our knowledge, these are the first such simulations in which multicomponent diffusion and interactions between dissociation and ionization of different species are consistently accounted for. Turbulence effects were represented by a standard k-epsilon model, both with and without an axisymmetric jet correction term and for several different choices of the turbulent Prandtl and Schmidt numbers Pr(t) and Sc(t). Simulations were performed for one FSH experiment and two BP experiments at different values of torch power P and argon flow rate W. The inflow profiles in the FSH simulations were adjusted to match P, W, and the experimental data slightly downstream of the torch exit as closely as possible. The same profile shapes were then used to match P and W for the BP simulations, for which data near the torch exit were not available. Swirl was neglected except in one of the FSH calculations, where it was found to have negligible effect, as expected. Best results were obtained with the axisymmetric jet correction term omitted and with Pr(t) = Sc(t) = 0.7. Agreement with the experimental data was then fair overall, but still showed systematic deviations and cannot be regarded as fully satisfactory. Possible reasons for the discrepancies are discussed.