The performance of the $k-{varepsilon}$ and $k-{omega}$ two-equation turbulence models was investigated in computational simulations of the neutrally stratified atmospheric boundary layer developing above various terrain types. This was achieved by using a proposed methodology that mimics the experimental setup in the boundary layer wind tunnel and accounts for a decrease in turbulence parameters with height, as observed in the atmosphere. An important feature of this approach is pressure regulation along the computational domain that is additionally supported by the nearly constant turbulent kinetic energy to Reynolds shear stress ratio at all heights. In addition to the mean velocity and turbulent kinetic energy commonly simulated in previous relevant studies, this approach focuses on the appropriate prediction of Reynolds shear stress as well. The computational results agree very well with experimental results. In particular, the difference between the calculated and measured mean velocity, turbulent kinetic energy and Reynolds shear stress profiles is less than ${pm}10%$ in most parts of the computational domain.