Based on preliminary theoretical analysis and numerical experiment, it is found that land surface heterogeneity plays an important role in the models turbulent flux calculation. In nearly neutral atmosphere conditions, variation coefficient of sub-scale roughness length, cell-average roughness, and reference height are main factors affecting the calculation of grid turbulent fluxes. The first factor has a determinant role on calculation deviation. The relative error generated by roughness heterogeneity could be more than 40% in some cases in certain areas (e.g., in vegetation-climate transition belt). Selecting a specific reference height may improve the calculation of turbulent flux. In stable or unstable atmosphere conditions, with sensible heat flux as an example, analysis shows that the discrepancy is correlated to the sub-grid distributions of mean wind velocity, potential temperature gradient between land surface and reference levels, and atmosphere stability near surface layer caused by the heterogeneity of land surface roughness. The calculation of turbulent flux is the most sensitive to stability in the above three factors. The above analysis shows that it is necessary to make a further consideration for the calculation deviation of the turbulent fluxes brought from land surface heteroge-neity in the present numerical models.