This paper describes a series of three‐dimensional simulations of shallow inefficient convection in the outer layers of the Sun. The computational domain is a closed box containing the convection–radiation transition layer, located at the top of the solar convection zone. The most salient features of the simulations are that: (i) the position of the lower boundary can have a major effect on the characteristics of solar surface convection (thermal structure, kinetic energy and turbulent pressure); (ii) the width of the box has only a minor effect on the thermal structure, but a more significant effect on the dynamics (rms velocities); (iii) between the surface and a depth of 1 Mm, even though the density and pressure increase by an order of magnitude, the vertical correlation length of vertical velocity is always close to 600 km; (iv) in this region the vertical velocity cannot be scaled by the pressure or the density scaleheight; this casts doubt on the applicability of the mixing length theory, not only in the superadiabatic layer, but also in the adjacent underlying layers; (v) the final statistically steady state is not strictly dependent on the initial atmospheric stratification.

The general consensus is that in order to reproduce the observed solar p-mode oscillation frequencies, turbulence should be included in solar models. However, until now there has not been any well-tested efficient method to incorporate turbulence into solar modeling. We present here two methods to include turbulence in solar modeling within the framework of the mixing length theory, using the turbulent velocity obtained from numerical simulations of the highly superadiabatic layer of the sun at three stages of its evolution. The first approach is to include the turbulent pressure alone, and the second is to include both the turbulent pressure and the turbulent kinetic energy. The latter is achieved by introducing two variables: the turbulent kinetic energy per unit mass, and the effective ratio of specific heats due to the turbulent perturbation. These are treated as additions to the standard thermodynamic coordinates (e.g. pressure and temperature). We investigate the effects of both treatments of turbulence on the structure variables, the adiabatic sound speed, the structure of the highly superadiabatic layer, and the p-mode frequencies. We find that the second method reproduces the SAL structure obtained in 3D simulations, and produces a p-mode frequency correction an order of magnitude better than the first method.