Direct Numerical Simulation of a Paradigmatic Urban Heat Island

In recent decades, global population growth has accelerated, and various projections indicate that it will continue rising steadily in the coming years. As a result, urban areas are expected to expand significantly, given their high population density. Whitin this evolving urban context, cities are particularly vulnerable to meteorological hazards due to the concentration of people, assets, infrastructure, and capital stocks [1], which has led to growing concern over urban climate issues and contributed to the rapid expansion of this research field [2].

Among the most studied phenomena in urban climate is the Urban Heat Island (UHI). As the term implies, UHIs are urbanized areas that are warmer than their surroundings, with a wide range of consequences, from impacts on human health to ecological changes.

Urban systems are inherently complex, and the UHI effect is largely influenced by the artificialization of the surface [3]. Heat produced by human activities and absorbed from solar radiation during the day is mainly stored within urban canopy elements and subsequently released during night time, making cities warm even during the night hours.

For this reason, the UHI phenomenon can be classified into diurnal and nocturnal types and can also be categorized based on the geographical setting of the city, since the surrounding areas interact with the urban canopy.

From its early documentation in London [4], the UHI phenomenon has presented a variety of scientific challenges, as it differs across cities and is driven by multiple interacting factors. As such, it has attracted considerable attention from disciplines including climatology, environmental science, architecture, materials science, energy engineering, urban planning, and public health [5].

Several studies in literature have examined the UHI phenomenon in real urban environments such as Hannover (Germany) [6], Toulouse (France) [2], and city across South Asia, for a review see the study of Kotharkar et al. [7] as some examples.

These studies provide valuable insight into how UHI manifests in diverse geographic and climatic contexts.

One of the main drivers of UHI effect is the urban canopy, which lies within the Urban Canopy Layer (UCL), the atmospheric layer extending from the ground up to the height of the dominant urban elements.

Understanding the complex dynamic that occurs within UCL is essential for accurately representing cities in weather and climate models. Since the urban microclimate is shaped by a range of physical processes occurring within the canopy layer, reliable urban parametrizations are essential for improving the performance of weather forecasting models.

At the neighborhood scale ranging from few meters up to 1 kilometer, variations in building form and density significantly influence airflow, thereby affecting the transport of momentum, heat, and pollutants within the urban environment [8]. However, in mesoscale atmospheric models such as, weather research and forecasting WRF model, buildings are not explicitly resolved. Instead, urban areas are represented through simplified, parameterized surface characteristics.

Building morphology plays a key role in these parameterizations. Urban models often describe city structures using morphological parameters such as the plan area index, ratio between the total buildings area and the total plan area, the frontal area index, the ratio between the building frontal area to the total plan area, and the average building height. Many parameterizations for airflow above cities are based on these and similar surface descriptors [1, 9-11].

Although few studies have attempted to simulate small portions of real urban, but simplified, environments, such as RANS simulations over parts of central London [12] on in Zhao old town [13], as in Basel [14], or wind tunnel experiments replicating selected districts [15, 16], most numerical investigations still adopt idealized geometries.

These configurations typically consist of simplified arrangements of cubic or rectangular blocks, either aligned or staggered, and often include elongated blocks to represent urban canyons. Such geometric simplifications are commonly used in both CFD simulations and wind tunnel experiments.

For instance, Xie and Castro [17] investigated flow over staggered wall-mounted cube of different heights using numerical simulations, Ricciardelli and Polimeno [18] explored airflow behavior in wind tunnels; and Zaki et al. [19] examined vertical wind profile within a canopy layer. Additional numerical studies of urban canopy flow on staggered and non-staggered cuboids are reported by Castro et al. [20].

Idealized configurations remain essential for identifying general flow behaviors and isolating the effects of key morphological parameters in urban environments.

Our study aligns with this line of research, aiming to develop a general understanding of the UHI phenomenon from a fluid dynamic perspective. To achieve this, we adopt a bottom-up approach, which aims to identify and analyze the fundamental physical processes governing urban atmospheric flows. This method emphasizes the role of key parameters, rather than relying on detailed reconstructions of specific urban configurations, which would constrain the applicability of the results to a narrow range of scenarios.

Since our analysis focuses on parameters related solely to urbanization, in our approach, the urban canopy is modelled as a rough surface composed of simplified hexahedral buildings arranged in a circular domain of radius R and not aligned with the computational grid. The urban layout is characterized by five key non-dimensional parameters: The average coverage frequency $\lambda$, the average density coverage $\rho$, the average aspect ratio of the planform area $A R_\pi$, and the average aspect ratio of the building heights $A R_z$ along with ε, a coefficient describing the radial variation of the city shape. As it will be shown these parameters effectively capture the main statistical features of cities and enable a generalized modelling approach. To the authors’ knowledge, this study presents the first investigation of the UHI under purely convective conditions, in the absence of external wind, within this paradigmatic setup using Direct Numerical Simulation.

The paper is organized as follows. Section 2 describes the algorithm employed to generate the idealized urban pattern. Section 3 outlines the numerical method and the governing equations. The flow topology and main statistics are described in Section 4, while Section 5 concludes the study with final remarks.

The flow is investigated within a slightly modified Rayleigh–Bénard Convection (RBC) setup. The governing equations are formulated under the Boussinesq approximation, in which density is assumed to be constant except for when it directly causes buoyant forces [21]. They are expressed in the following non dimensional form:

$\frac{\partial u_i}{\partial x_i}=0$            (14)

$\frac{\partial u_i}{\partial t}+\frac{\partial u_i u_j}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\sqrt{\frac{P r}{R a}} \frac{\partial^2 u_i}{\partial x_j \partial x_j}+\vartheta \delta_{i_3}$              (15)

$\frac{\partial \vartheta}{\partial t}+\frac{\partial \vartheta u_j}{\partial x_j}=\frac{1}{\sqrt{\operatorname{RaPr}}} \frac{\partial^2 \vartheta}{\partial x_i \partial x_i}$               (16)

The variables $u_i, p$, and $\vartheta$ denote the velocity, pressure, and the temperature fields respectively. The indices i,j =1, 2, 3, where $u_1=u$, and $u_2=v$ denotes the horizontal velocity components, $u_3=w$ is the vertical component, and $\delta_{i, j}$ represents the Kronecker delta. The dimensionless form Eqs. (14)-(16) is obtained by using the cell height H , the temperature difference $\Delta \Theta$ between the maximum temperature within the city and the temperature of the top plate and the free-fall velocity $U_{f f}=\sqrt{g \beta \Delta \Theta H}$ as characteristic scales. The Prandtl number is given by the ratio between the kinematic viscosity $v$ and thermal diffusivity $\alpha$, i.e., $\operatorname{Pr}= v / \alpha$ while the Rayleigh number is given by $R a= g \beta \Delta \Theta H^3 / v \alpha$ where β is the thermal expansion coefficient and g is the gravitational acceleration.

3.2 Numerical method

The Direct Numerical Simulation of this idealized city pattern is performed using the open-source code CaNS [22]. Spatial discretization is based on a second order finite difference scheme on a staggered grid, while time integration is achieved by a three-step Runge-Kutta scheme with a Courant–Friedrichs–Lewy (CFL) condition of 0.6.

The simulation employs an immersed boundary method (IBM), specifically the formulation by Yang and Balaras [23], to handle the complex urban geometry. To this end, a signed distance field $\phi$, (SDF) of the created geometry is provided as input to the solver with values defined at each grid point in the domain. This scalar field enables a representation of the solid boundaries: positive values of $\phi$ indicate fluid regions, negative values correspond to solid regions and the fluid–solid interface is defined by the iso-surface $\phi=0$.

This representation is essential for identifying the forcing points, grid nodes closest to the immersed surface, as well as the interpolation points, used to reconstruct the velocity field at those locations. In IBM approaches, an additional force term is added into the right-hand side of the Navier-Stokes equation to enforce the desired velocity boundary condition at the immersed surface. In our case, the IBM consists of four main steps. The first two are executed only once at the beginning of the simulation, since the boundaries are stationary. In this initialization phase, the solver determines the forcing points, defined above. Then, for each forcing point, it determines three interpolation points, which are later used to reconstruct the velocity field at those locations.

The remaining steps are performed at each time step. Hence, the velocity field is interpolated at the previously identified interpolation points, and the corresponding forcing term is computed and applied to the forcing points to enforce the desired boundary condition at the fluid-solid interface.

This simulation is performed for a fixed Prandtl number Pr = 0.7, corresponding to air, and for Rayleigh number, $R a=10^8$. The computational domain is a rectangular box of size $(L x, L y, L z)=(8,8,1.05)$, where x, y are the horizontal axes and z is the wall normal vertical coordinate. It is worth noting that the reference length used for non-dimensionalization is H = 1, corresponding to vertical distance between the ground surface, at the base of the buildings and the top plate.

The domain is discretized using $(n x, n y, n z)=(1600,1600,263)$ set grid points. In the vertical direction, a hyperbolic tangent stretching is applied to ensure the correct resolution near the walls. 

Boundary conditions for the velocity field are defined as follows: periodic boundary conditions are applied at the lateral directions, while no-slip boundary condition is used at the bottom plate and free-slip boundary condition is set at the top plate. A schematic representation of the domain can be seen in Figure 3.

Regarding the temperature field, an isothermal condition is imposed at the top plate, with $\vartheta_{\text {cold }}=-0.5$, while lateral sides are treated as periodic. At the bottom plate, a non-uniform temperature distribution is applied: within the urban pattern as shown in Figure 4, a Gaussian profile with a peak value of $\vartheta_{\text {hot }}=0.5$ is prescribed, whereas a constant temperature is imposed outside this region. The maximum temperature difference across the fluid layer is thus defined as $\Delta=\vartheta_{\text {hot }}- \vartheta_{\text {bottom }}=1$.

$T_{\text {top }}=T_{\text {atmosphere }} \leq T_{\text {country }} \leq T_{\text {city }}$

Figure 3. Schematic representation of the computational domain and boundary conditions

Figure 4. Skyline of the created city pattern

The simulation parameters are summarized in Table 2.

After an initial transient, velocity, pressure and temperature fields are collected when the system reaches a statistically steady state characterized by an average Nusselt number of Nu = 15 computed at the top plate, and with a time interval larger than the time eddy turnover given by $T_{ {free fall }}=H / U_{f f}$. 

Table 2. Simulation parameters

Notes: $\Delta \tau$ is the time interval between consecutive fields used for statistical averaging; T is the total simulation time

Recalling the final set of parameters used for the city generation, as reported in Table 1, we can identify a non-dimensional roughness thickness $k_z^*$,  given by the ratio between the thermal boundary layer thickness $\delta_t=H /(2 N u)$, and the mean buildings height, that provides a measure of how strongly the urban geometry interacts with the thermal boundary layer. In our case, we find $k_z^* \sim 0.88$.

UHI is one of the most widely studied phenomena in urban climatology over the past decades. In this study, we have developed a general framework to represent a variety of idealized urban patterns using only five parameters. These parameters capture key morphological characteristics, such as coverage density, mean height, and spatial distribution, that are measurable from real urban layouts.

For the first time, a Direct Numerical Simulation of such an idealized configuration has been performed, consisting of 176 buildings, with the aim of exploring the fundamental mechanisms governing heat transfer under buoyancy-driven flow conditions. The proposed setup provides a highly controlled and repeatable environment to isolate the influence of urban morphology on thermal processes.

While many other studies focus on specific neighborhoods, the present work adopts an idealized and general approach for an entire city pattern. This generalization allows for a more fundamental understanding of the physical mechanism, free from the geometric and contextual constraints of a specific urban case.

Findings obtained through this approach can be generalized and subsequently applied to more realistic and site-specific configurations

Future developments of this work may include the introduction of a mean wind to assess the combined effect of shear and buoyancy driven flows. Additionally, the influence of both physical parameters, such as $k_z^*$ and morphological features of the city could be systematically varied to better understand the dominant mechanisms that control thermal exchanges and to identify effective strategies for urban climate mitigation. Moreover, given the high computational cost of DNS, future investigations could benefit from the use of Large Eddy Simulations (LES) to explore a wider range of configurations within reasonable time frames and, above all, reduce the time costs.

In summary, this study provides a novel and flexible approach that can serve as a reference for future numerical investigations on UHIs.