A 3-D numerical method is developed to investigate the spatial distribution of surface fluxes over heterogeneous surfaces in (semi-)arid regions. Quantifying the effects of changes in the momentum, thermal and moisture roughness lengths on the airflow and fluxes in the ABL is important for water resources management and local climate studies. The governing equations and turbulence models are modified to include the effects of atmospheric stability conditions on the airflow. The turbulent airflow in ABL is simulated based on the Unsteady Reynolds-Averaged Navier-Stokes (URANS) approach to understand the air flow over the non-homogeneous surfaces from dry land through the water surface and vice versa. The model can be used to study airflow in neutral and non-neutral ABL over complex and non-homogeneous surfaces. The model results were used to investigate the flow parameters and (heat) flux variations over small water surfaces considering its surrounding conditions.