We present the first numerical simulations of the solar interior to exhibit a tachocline consistent with the Gough and McIntyre (1998) model. We find nonlinear, axisymmetric, steady-state numerical solutions in which: (1) a large-scale primordial field is confined within the radiation zone by downwelling meridional flows that are gyroscopically pumped in the convection zone (2) the radiation zone is in almost-uniform rotation, with a rotation rate consistent with observations (3) the bulk of the tachocline is magnetic free, in thermal-wind balance and in thermal equilibrium and (4) the interaction between the field and the flows takes place within a very thin magnetic boundary layer, the tachopause, located at the bottom of the tachocline. We show that the thickness of the tachocline scales with the amplitude of the meridional flows exactly as predicted by Gough and McIntyre. We also determine the parameter conditions under which such solutions can be obtained, and provide a simple explanation for the failure of previous numerical attempts at reproducing the Gough and McIntyre model. Finally, we discuss the implications of our findings for future numerical models of the solar interior, and for future observations of the Sun and other stars.

The solar convection zone exhibits a strong level of differential rotation, whereby the rotation period of the polar regions is about 25-30% longer than the equatorial regions. The Coriolis force associated with these zonal flows perpetually "pumps" the convection zone fluid, and maintains a quasi-steady circulation, poleward near the surface. What is the influence of this meridional circulation on the underlying radiative zone, and in particular, does it provide a significant source of mixing between the two regions? In Paper I, we began to study this question by assuming a fixed meridional flow pattern in the convection zone and calculating its penetration depth into the radiative zone. We found that the amount of mixing caused depends very sensitively on the assumed flow structure near the radiative--convective interface. We continue this study here by including a simple model for the convection zone "pump", and calculating in a self-consistent manner the meridional flows generated in the whole Sun. We find that the global circulation timescale depends in a crucial way on two factors: the overall stratification of the radiative zone as measured by the Rossby number times the square root of the Prandtl number, and, for weakly stratified systems, the presence or absence of stresses within the radiative zone capable of breaking the Taylor-Proudman constraint. We conclude by discussing the consequences of our findings for the solar interior and argue that a potentially important mechanism for mixing in Main Sequence stars has so far been neglected.