Consistency and stability are two essential ingredients in the design of numerical algorithms for partial differential equations. Robust algorithms can be developed by incorporating nonlinear physical stability principles in their design, such as the entropy production inequality (or second law of thermodynamics), rather than by simply adding artificial viscosity. In the context of two-equation turbulence models we introduce space-time averaged variables, the essential concept which enables identification of an appropriate set of conservation variables. From these, the correct concept of generalized entropy and a set of entropy variables can be defined which leads to a symmetric system of advective diffusive equations, as established by Mock and Godunov. Positivity and symmetry of the equations require certain constraints on the turbulence diffusivity coefficients and the turbulence source terms. With these, we are able to design entropy producing formulations of two-equation turbulence models and, in particular, the k-epsilon model, and numerical formulations that inherit these properties. The accuracy of the original k-epsilon model is maintained and we automatically gain computational stability and robustness due to the guaranteed entropy production property.