Abstract Multiple zonal jets are investigated with a two-level primitive equation model on the sphere in which both baroclinicity and planetary radius are varied. As in the case for a two-layer quasigeostrophic model on a β-plane channel, it is found both that the Rhines scale successfully predicts the meridional scale of the multiple zonal jets, and that these jets are maintained in part by an eddy momentum flux divergence associated with slow baroclinic waves at the interjet minimum. A scaling analysis suggests that njets∝ (a/θm)1/2, with the constraints ζe ≡ 8 sin2f (θm/▵θ ) > 1 and njets ≥ 1, where njets is the number of the jets, a the planetary radius, θm one-half of the pole-to-equator potential temperature difference, ξe the supercriticality of the two-layer Phillips model, Δθ the potential temperature difference between the two levels, and ϕ the latitude. The number of jets simulated by the model agrees with this scaling, provided that Ljet ≤ a, where Ljet is the jet scale. In model runs with a large planet where multiple zonal jets exist, the time–mean eddy heat flux is found to be consistent with the diffusive picture of Held and Larichev. In contrast, for the model runs with the planetary size equal to that of Earth, baroclinic adjustment is found to be more relevant. These results are consistent with the finding that in the large-planet (Earth-like) model runs, the jet/eddy scale is smaller than (comparable to) the corresponding planetary radius.