The solar convection zone exhibits a differential rotation with radius and latitude that poses major theoretical challenges. Helioseismology has revealed that a smoothly varying pattern of decreasing angular velocity Ω with latitude long evident at the surface largely prints through much of the convection zone, encountering a region of strong shear called the tachocline at its base, below which the radiative interior is nearly in uniform solid body rotation. Helioseismic observations with MDI on SOHO and with GONG have also led to the detection of significant variations in Ω with 1.3 yr period in the vicinity of the tachocline. There is another shearing layer just below the solar surface, and that region exhibits propagating bands of zonal flow. Such rich dynamical behavior requires theoretical explanations, some of which are beginning to emerge from detailed 3-D simulations of turbulent convection in rotating spherical shells. We discuss some of the properties exhibited by such numerical models. Although these simulations are highly simplified representations of much of the complex physics occurring within the convection zone, they are providing a very promising path for understanding the solar differential rotation and its temporal variations.