When wind turbines are arranged in large wind farms, their efficiency decreases significantly due to wake effects and to complex turbulence interactions with the atmospheric boundary layer (ABL) [1]. For large wind farms whose length exceeds the ABL height by over an order of magnitude, a “fully developed” flow regime may be established [1, 2, 3]. In this asymptotic regime, the changes in the stream-wise direction are small compared to the more relevant vertical exchange mechanisms. Such a fully developed wind-turbine array boundary layer (WTABL) has recently been studied [2] using Large Eddy Simulations (LES) under neutral stability conditions. The simulations showed the existence of two log-laws, one above (characterized by: uhi ∗ , zhi o ) and one below (ulo ∗ , zlo o ) the wind turbine region. This enabled the development of more accurate parameterizations of the effective roughness scale for a wind farm. Now, a suite of Large Eddy Simulations, in which wind turbines are modeled as in [2] using the classical drag disk concept are performed, again in neutral conditions but also considering temperature. Figure 1 shows a schematic of the geometry of the simulation. The aim is to study the effects of different thermal ABL stratifications, and thus to study the efficiency and characteristics of large wind farms and the associated land-atmosphere interactions for realistic atmospheric flow regimes. Such studies help to unravel the physics involved in extensive aggregations of wind turbines, allowing us to design better wind farm arrangements. As a first step, temperature is treated in a passive mode, allowing us to focus the study on the influence of a large WFABL into the scalar fluxes. By considering various turbine loading factors, surface roughness values and different atmospheric stratifications, it is possible to analyze the influence of these parameters on the induced surface roughness, and the sensible heat roughness length. These last two parameters can be used to model wind turbine arrays in simulations of atmospheric dynamics at larger (regional and global) scales [4], where the coarse meshes used do not allow to account for the specifics of each wind turbine. Results from different sets of simulations are presented, for which also the corresponding effective roughness length-scales can be determined. The results also help our understanding of how wind turbines affect scalar transport processes in the turbine wakes. By using a simple drag disk approach for modeling the wind turbines, it is found that the surface heat flux inside the thermal wind-turbine array boundary layer is increased. This is the result of two competing effects: (1) a major increase on u∗,hi; (2) a smaller decrease due to lower u∗,lo near the ground.