Waves breaking in the shallow surf zone near the shoreline inject turbulence into the water column that may reach the bed to suspend sediment. Breaking-wave turbulence in the surf zone is, however, poorly understood, which is one of the reasons why many process-based coastal-evolution models predict coastal change during severe storms inadequately. Here, we use data collected in two natural surf zones to derive a new parameterization for the stability function C-mu that determines the magnitude of the eddy viscosity nu(t) in two-equation turbulent-viscosity models, nu(t) = C(mu)k(2)/epsilon, where k is turbulent kinetic energy and epsilon is the turbulence dissipation rate. In both data sets, the ratio of turbulence production to dissipation is small (approximate to 0.15), while vertical turbulence diffusion is significant. This differs from assumptions underlying existing C-mu parameterizations, which we show to severely overpredict observed C-mu for most conditions. Additionally, we rewrite our new C-mu parameterization into a formulation that accurately reproduces our Reynolds-stress based estimates of turbulence production. This formulation is linear with strain, consistent with earlier theoritical work for large strain rates. Also, it does not depend on epsilon and can, therefore, also be applied in one-equation turbulent-viscosity models. We anticipate our work to improve turbulence modeling in natural surf zones and to eventually lead to more reliable predictions of coastal evolution in response to severe storms.