scale heights counted from the upper boundary. The very large superadiabatic gradients near the upper boundary are mainly caused by the solid boundary conditions. In the middle region of convection zone, the small superadiabatic gradient could be affected by the numerical errors. (see text for details)[]{data-label="sags"}](./figures/sag-ms.eps "fig:"){width="8.2cm"}

After a long period of thermal relaxation, the flow reaches a statistically steady state and its average properties become stationary. One criterion of thermal relaxation is that the input energy flux from the bottom is balanced by the outgoing energy flux through the top surface. In the present calculations, the deviations of the averaged vertical total flux from the input energy flux are less than $2$% everywhere (except Case I1 which has a bound of $6$%). The maximum discrepancies occur near the unstable-stable interfaces (where the cross-correlation of temperature and vertical velocity changes its sign).

In expressing the results, the overline denotes a combination of horizontal and temporal averaging, the prime ($'$) denotes the deviation from the mean, and the double prime ($''$) denotes the rms fluctuation from the mean. In some of the figures, the locations of integral pressure scale heights counted from the upper boundary are also shown. For example, in 