  Renzini (1987) wrote an influential critique of mixing-length theory (MLT) as
used in stellar evolution codes, and concluded that three-dimensional (3D)
fluid dynamical simulations were needed to clarify several important issues. We
have critically explored the limitations of the numerical methods and conclude
that they are approaching the required accuracy. Implicit large eddy
simulations (ILES) automatically connect large scale turbulence to a Kolmogorov
cascade below the grid scale, allowing turbulent boundary layers to remove
singularities that appear in the theory. Interactions between coherent
structures give multi-modal behavior, driving intermittency and fluctuations.
Reynolds averaging (RA) allows us to abstract the essential features of this
dynamical behavior of boundaries which are appropriate to stellar evolution,
and consider how they relate static boundary conditions (Richardson,
Schwarzschild or Ledoux). We clarify several questions concerning when and why
MLT works, and does not work, using both analytical theory and 3D high
resolution numerical simulations. The composition gradients and boundary layer
structure which are produced by our simulations suggest a self-consistent
approach to boundary layers, removing the need for ad hoc procedures for
'convective overshooting' and `semi-convection'. In a companion paper we
quantify the adequacy of our numerical resolution, determine of the length
scale of dissipation (the `mixing length') without astronomical calibration,
quantify agreement with the four-fifths law of Kolmogorov for weak
stratification, and extend MLT to deal with strong stratification.
