The integrated hydrodynamic thermal-capillary model (IHTCM) of Czochralski growth for large-diameter silicon crystals is extended to include a k-epsilon model for turbulence in the melt implemented in a form appropriate for capturing the transition to nearly laminar flow near solid boundaries. Calculations are presented for buoyancy-driven flow alone and for the flow driven by a combination of crystal and crucible rotation, buoyancy and surface tension gradients. These results predict the enhancement in the heat and mass transfer seen in experiments with increased crucible rotation rate, which is not predicted by laminar flow simulatons. The computed temperature fields and interface shapes compare well with measurements reported before (Kinney, Bornside, Brown and Kim, J. Crystal Growth 126 (1992) 413). The use of the k-epsilon/IHTCM for optimization of operating conditions is demonstrated by calculations for varying crystal and crucible rotation rates using an objective function that attempts to optimize oxygen concentration in the crystal, to minimize the radial variation of oxygen and to reduce the magnitude of the thermoelastic stress.