d}$ is defined in (\[fd\]) and the total flux is $F_{\rm t}=F_{\rm e}+F_{\rm
k}+F_{\rm v}+F_{\rm d}$. A more detailed discussion of these fluxes is given in Section \[sec:flx\].

As an important indicator of convective instability, the superadiabatic gradient $\Delta\nabla$ is given in Fig. \[sags\] for all the cases. Here we focus on cases I1–I3. Highly superadiabatic regions occur near the upper stable-unstable interfaces, the systems are substantially unstable. The sharp jumps may be caused by the lack of enough room for upward motions. The same phenomenon was found in [@singh98]’s work (see Fig.5 therein). In the convection zone, $\Delta\nabla$ is not as small as that of the 1D reference model. It is caused by the enhanced total flux. The slightly negative $\Delta\nabla$ in the lower part of the convection zone ($0.4 < z < 0.6$) indicates $\Delta\nabla$ is sensitive to even small numerical inaccuracy. A higher-order scheme with denser mesh may eliminate such problem.

Figure \[vfs\] shows examples of the velocity fields projected onto a vertical plane. We can see from these figures that the turbulent flows are dominated by large eddies whose sizes are comparable to at least a press scale height. The small scale flows are mostly associated 