Introduction

Should engineers and urban planners care about CFD in cities?

Heatwaves are no longer rare events. With three major heatwaves in may, june and july 2026 in Europe, they are becoming the defining climate challenge of our era, turning dense urban environments into heat traps where temperatures can exceed surrounding rural areas by 5 to 10°C. This phenomenon, the Urban Heat Island (UHI) effect, kills thousands of people every summer and strains energy grids to their breaking points.

The question is: can we simulate it accurately enough to design cities with a better resilience to heatwaves ?

A 2022 study published in Sustainability by Y. Qu, M. Milliez, L. Musson-Genon, and B. Carissimo answers with a resounding yes using code_saturne, the open-source CFD solver developed by EDF R&D. Their work simulates the full diurnal thermal cycle of a real neighborhood in the city center of Toulouse, France, validates it against an exceptional experimental dataset, and then uses it to reveal how urban heat fundamentally reshapes pollutant dispersion.

In this article, we will walk you through their methodology, their model and most importantly, their experimental campaigns and code_saturne results.

The core problem: why neutral atmosphere assumptions are dangerous

Most micro-meteorological studies of urban flow and pollution dispersion assume a neutral atmosphere. This means the temperature variations are ignored and buoyancy forces are set to zero. This is computationally convenient but physically wrong, especially during summer days when sun-baked rooftops reach temperatures far above ambient air. Qu et al. demonstrate that this assumption is not just an approximation; it is an error that:

• Misrepresents surface temperatures by tens of degrees
• Underestimates turbulent kinetic energy near rooftop level
• Completely alters predicted pollutant plume shapes

The solution they propose is a fully coupled thermo-radiative-convective model, implemented in code_saturne's atmospheric module, applied to a real district of Toulouse during the CAPITOUL (Canopy and Aerosol Particles Interactions in TOulouse Urban Layer) field campaign.

Presentation of the resolved physics by the model

The physics of urban heat

In cities, heat transfer mainly involves three phenomena that cannot be treated independently:

• Radiative heat transfer: solar radiation heating surfaces, surfaces re-emitting thermal infrared radiation, buildings exchanging radiation with each other and the sky.
• Convective heat transfer: wind carrying heat through street canyons and around buildings, with the critical heat transfer coefficient varying spatially at every wall patch.
• Conductive heat transfer: heat diffusing through walls, rooftops, and the ground subsurface.

The study central contribution is encoding all three into a single, consistent simulation framework.

The surface energy balance

For every surface in the urban domain, every wall patch, every roof tile, every road segment, the model solves a Surface Energy Balance (SEB):

$$Q_{cond} + Q_H = Q^*$$

where:

• $Q_{cond} (W.m^{-2})$ is the conductive heat flux linking the surface temperature to the internal building temperature or deep soil temperature.
• $Q_H (W.m^{-2})$ is the sensible heat flux exchanged with the local airflow.
• $Q^* (W.m^{-2})$ is the net radiative flux (shortwaves from the sun + longwave from the buildings thermal radiations).

The net radiative flux $Q^*$ is itself decomposed into its shortwave and longwave components:

$$Q^* = S^* + L^*$$

where:

• $S^* = (1 - \alpha)(S_D + S_f + S_e)$ $(W.m^{-2})$ is the net shortwave radiative flux: the solar energy actually absorbed by the surface after reflection, with $\alpha$ the surface albedo and $S_D$, $S_f$, $S_e$ the direct, atmosphere-diffuse, and environment-reflected solar components respectively.
• $L^* = \varepsilon(L_a + L_e - \sigma T_{sfc}^4)$ $(W.m^{-2})$ is the net longwave radiative flux: the balance between the longwave radiation received from the sky ($L_a$) and surrounding buildings ($L_e$), and the longwave radiation emitted by the surface itself as a gray body ($\varepsilon \sigma T_{sfc}^4$).

Note that anthropogenic heat flux (heat emitted by human activity) and latent heat flux are neglected in this study. This simplification would eventually need to be lifted for a fully general model of a living city, but is appropriate for the dry summer conditions of the CAPITOUL campaign.

The sun emits short-wave radiations, absorbed and re-emitted as long waves radiations by building: these are encoded in $Q^*$. The sensible heat $Q_H$ is the heat exchanged with local airflow, that will then advect (transport) it locally. Finally, the buildings and soil have their own temperature that interacts with air by conductive processes $Q_{cond}$. Picture taken from [1]. Anthropogenic heat and latent heat (heat released or absorbed due to phase transitions) are neglected here.

This equation is not solved once for an average surface. It is solved at every surface cell in the mesh, at every timestep, with $Q_H$ computed from the local 3D wind field and $Q^*$ computed from the 3D radiative transfer solver. This is what makes the model genuinely three-dimensional.

From heat balance to surface temperature $T_{sfc}$

One of the most elegant aspects of this paper is the hybrid thermal scheme: ground temperature and building surface temperature are treated with different models, each chosen to best match the underlying physics.

Ground Temperature: The Force-Restore Method

For the soil, the authors use the classical force-restore method (Deardorff 1978), widely used in meteorological models:

$$\frac{\partial T_{sfc}}{\partial t} = \sqrt{2\omega} \, \mu_t (L^* + S^* - Q_H) - \omega(T_{sfc} - T_{g/b})$$

where:

• $T_{sfc}$ is the surface temperature on the ground.
• $\omega=2\pi/86400$ ($s^{-1}$).
• $L^* (W.m^{-2})$ is the net long-wave radiative flux introduced above.
• $S^* (W.m^{-2})$ is the net short-wave radiative flux introduced above.
• $Q_H (W.m^{-2})$ is the sensible heat flux.
• $\mu_t (J.m^{-2}.s^{-0.5}.K^{-1})$ is the thermal admittance of the soil.
• $T_{g/b}(K)$ is the deep soil (or internal building) temperature acting as a restoring temperature.

The "force" term (the first term on the right) drives the surface temperature in response to the net energy forcing. The "restore" term (the second) prevents the surface from drifting too far from its deep equilibrium temperature on diurnal timescales. This is an elegant, physically motivated, and computationally cheap way to handle soil thermal inertia.

Building Surface Temperature: The Wall Thermal Model

For building walls and rooftops, a single-layer conduction model is used. The energy balance at a wall surface is:

$$\frac{\lambda}{e}(T_{sfc} - T_{int}) + h_f(T_{sfc} - T_a) = L^* + S^* = \varepsilon(L_a + L_e - \sigma T_{sfc}^4) + (1 - \alpha)(S_D + S_f + S_e)$$

Let us decode this equation term by term; it is the heart of the thermal model:

Left-hand side:

• $\frac{\lambda}{e} \left( T_{sfc} - T_{int} \right)$: conductive heat flux through the wall, where $\lambda (W. K^{-1}.m^{-1})$ is the average thermal conductivity and $e (m)$ is the wall thickness. This links the surface temperature to the internal building temperature $T_{int}$.
• $h_f(T_{sfc}- T_a)$ is the convective heat flux to the external air, where $h_f (W.m^{-2}.K^{-1})$ is the local heat transfer coefficient and $T_a$ is the external air temperature.

Right-hand side:

• $L^*=\epsilon (L_a + L_e - \sigma T_{sfc}^4$) : net long-wave radiative flux. $L_a$ is the atmospheric downwelling long-wave radiation, $L_e$ is the long-wave radiation from surrounding building surfaces (the crucial wall-to-wall radiation), and $\epsilon \sigma T_{sfc}^4$ is the surface's own emission as a gray body.

• $S^*=(1-\alpha)(S_D+S_f+S_e)$: absorbed short-wave radiation, where $\alpha$ is the surface albedo and $S_D$, $S_f$, $S_e$ are the direct solar, atmosphere-diffuse, and environment-reflected solar components respectively.

This equation is solved implicitly for $T_{sfc}$ at each surface patch and each timestep, coupling the radiation solver, the CFD flow solver, and the wall conduction model simultaneously. Now, we need to explain how the radiative transfer terms $L_a$, $L_e$, $S_D$, $S_f$, and $S_e$ are computed.

Radiative transfer $Q^*$

Solving the full 3D problem with the Discrete Ordinates Method (DOM)

The radiative fluxes are computed using the Discrete Ordinates Method (DOM), which solves the Radiative Transfer Equation (RTE) for gray, non-diffusive, semi-transparent media. The key insight of the DOM implementation in Code_Saturne is that the RTE is solved in the entire fluid domain, not just at solid faces as would be the case with view factor methods. This means:

• The radiation field is computed volumetrically, making the approach naturally extensible to non-transparent media (fog, aerosol layers, pollution). It is however assumed here that the atmosphere is fully transparent to such temperature radiations.
• The spatial discretization uses the same mesh as the CFD solver: no interpolation between different grids.
• The angular discretization uses 32 or 128 discrete directions, depending on the required accuracy.

Both short-wave (solar) and long-wave (thermal infrared) radiation are treated separately, which is physically essential: building surfaces absorb solar radiation with a spectral albedo $\alpha$ and emit thermal radiation with emissivity $\epsilon$, and these two quantities are independent material properties.

Shortwave radiation: how the sun heats the city

The net shortwave flux absorbed by a surface is then:

$$S^* = (1 - \alpha)(S_D + S_f + S_e)$$

A key modeling challenge is the computation of shadows: $S_D$ is zero on a surface that is shaded by a neighboring building. The DOM solver, operating on the full 3D mesh, naturally handles this by propagating the solar beam through the domain and letting the geometry block it. This is one of the reasons why the paper uses the same mesh for both the CFD and radiation calculations: shadow boundaries are resolved at the same spatial resolution as the flow field (down to 0.5 m near the center of the domain).

In Toulouse, the authors further classified building surfaces into four albedo categories based on wall paint color, observed from Google Maps imagery: rose, gray, whitewash, and white. This reflects a physical reality: a white-painted wall in direct sunlight absorbs dramatically less solar energy than a dark brick wall and demonstrates that even a seemingly minor modelling detail (wall color) has a measurable effect on simulated surface temperatures.

Model visualization of $S^\uparrow$, to be compared and calibrated with experimental data.

The upward solar flux leaving a surface, measurable by a pyranometer on the mast, is:

$$S^\uparrow = \alpha (S_D + S_f + S_e)$$

This quantity is directly compared to measurements in the CAPITOUL validation, with very good agreement ($\sim 20\ W.m^{-2}$ overestimation at noon).

Longwave radiation: buildings radiating to each other

In parallel, every surface emits thermal (longwave, infrared) radiation proportional to the fourth power of its temperature. The net longwave flux at a surface is:

$$L^* = \varepsilon \left( L_a + L_e - \sigma T_{sfc}^4 \right)$$

where $L_a$ is the downwelling longwave radiation from the sky (treated as a blackbody at effective sky temperature), $L_e$ is the longwave radiation received from surrounding building surfaces, and $\varepsilon \sigma T_{sfc}^4$ is the surface's own gray body emission. The wall-to-wall longwave exchange term $L_e$ is what makes the DOM approach so powerful: it automatically accounts for the fact that a sun-baked south-facing wall at 60°C continuously irradiates the opposite north-facing wall across the street, keeping it warmer than it would be in isolation.

The outgoing longwave flux leaving a surface upward is:

$$L^\uparrow = \varepsilon \sigma T_{sfc}^4$$

This makes it extremely sensitive to surface temperature: a 3°C error in $T_{sfc}$ at 320 K propagates to approximately $20\ W.m^{-2}$ error in $L^\uparrow$. The simulation shows an underestimation of about $20\ W.m^{-2}$ during the night and an overestimation of up to $100\ W.m^{-2}$ at noon, both directly traceable to the corresponding surface temperature biases.

Brightness temperature: what infrared cameras actually measure

A key innovation in the validation strategy is the use of brightness temperature $T_{br}$ rather than actual surface temperature. Airborne infrared cameras measure radiance (not temperature directly) and it depends on both the true surface temperature and the emissivity. For a gray body surface in an urban canyon, accounting for infinite reflections between canyon walls, the brightness temperature is:

$$T_{br} = \left( \varepsilon T_{sfc}^4 + \frac{(1-\varepsilon) L^\downarrow}{\sigma} \right)^{1/4}$$

where $\varepsilon$ is the long-wave emissivity, $\sigma = 5.67 \times 10^{-8}$ $(W.m^{-2}.K^{-4})$ is the Stefan-Boltzmann constant, and $L^\downarrow$ $(W.m^{-2})$ is the incident long-wave radiation. The brightness temperature is always lower than the true surface temperature when $\varepsilon < 1$. The paper uses a commonly accepted approximation that holds well when the reflected longwave term $(1-\varepsilon)L^\downarrow/\sigma$ is small compared to $\varepsilon T_{sfc}^4$ i.e. when $T_{sfc}$ is large or $(1−\varepsilon)$ is small:

$$T_{br} \approx T_{sfc} \cdot \varepsilon^{1/4}$$

The paper shows that on a summer afternoon, this difference reaches approximately 3°C on rooftops: not negligible when validating against airborne infrared images. Comparing raw surface temperatures to IRT measurements without this correction would systematically overestimate temperatures.

The diurnal cycle: getting temperature right through the day

For the full diurnal simulation of 15 July 2004, the model reproduces the observed brightness temperatures at fixed IRT sensor positions around the Alsace-Lorraine and Pomme streets. Key findings:

• Street surfaces: the force-restore method for the ground works well; the main source of error is shadow timing: the model's approximated geometry sometimes misses the exact moment of sun/shade transitions.
• West-facing walls: slight overestimation in the morning (0600–1200 UTC), which the paper attributes to possible sensor underestimation. A similar finding was reported by Hénon et al. (2012) using the SOLENE model.
• Rooftops: excellent agreement before sunrise and from late afternoon onward; a morning advance in warming due to non-modeled small obstacles (chimneys, ventilation stacks) that would shade the sensor's detected surface during sunrise.

This last point is important: small unmodeled architectural details matter. A chimney that casts a shadow on a sensor cannot be ignored when validating at the single-cell level.

Conductive heat transfer $Q_{cond}$: what lies beneath the surface

Why conduction cannot be ignored

In the surface energy balance $Q_{cond} + Q_H = Q^*$, the conduction term $Q_{cond}$ represents the heat flux that flows into or out of the solid material (through a wall, a rooftop, or the ground) rather than into the air. It plays three physically critical roles:

• Thermal inertia: a massive concrete wall stores the solar energy it absorbs during the day and releases it slowly at night, keeping urban surfaces warm long after sunset. This is a major contributor to the nocturnal urban heat island.
• Coupling between exterior and interior: the exterior surface temperature $T_{sfc}$ is not independent of what happens inside the building. A well-insulated building will have a very different exterior surface temperature profile than an uninsulated one, even under identical radiative forcing.
• Stabilizing the energy balance: without a conduction sink, all absorbed radiative energy would be forced into the air by convection, leading to wildly overestimated surface temperatures (as the "radiation only" simulation clearly demonstrates).

Conduction through walls: the single-layer model

For building walls and rooftops, the model uses a single-layer approximation of the Fourier heat conduction equation. In steady-state through a wall of thickness $e$ and thermal conductivity $\lambda$, the conductive heat flux is:

$$Q_{cond} = \frac{\lambda}{e}(T_{sfc} - T_{int})$$

This is the linearized form of Fourier's law $Q_{cond} = -\lambda \nabla T$, applied across the wall thickness. The wall thermal resistance is $R = e/\lambda$ $(m^2.K.W^{-1})$, and the conductive flux is proportional to the temperature difference between the exterior surface $T_{sfc}$ and the interior building temperature $T_{int}$.

The thermal parameters used in the Toulouse simulation are derived from [2] and are summarized in the paper: for the brick and stone walls typical of the Toulouse city center, representative values are $\lambda \approx 1.0\ W.K^{-1}.m^{-1}$ and $e \approx 0.5\ m$, giving a thermal resistance of $R \approx 0.5\ m^2.K.W^{-1}$.

Tracking internal building temperature $T_{int}$

The internal building temperature $T_{int}$ is not a fixed boundary condition: it evolves during the diurnal cycle. The paper uses an incremental-adjustment method:

$$T_{int}^{n+1} = T_{int}^{n-1}\left(\frac{\tau - \Delta t}{\tau}\right) + T\left(\frac{\Delta t}{\tau}\right)$$

where $\tau = 86400\ s$ is the number of seconds in a day, $\Delta t$ is the timestep, and $T$ is the average outside building surface temperature computed at the current timestep $n$. This equation describes a first-order relaxation: the interior temperature tracks the exterior surface temperature but with a lag of order $\tau$ (one full day). Physically, this represents the slow thermal response of the building mass: the interior warms up gradually through the day as heat conducts inward through the walls, and cools slowly through the night.

For initialization, the internal temperature can be set equal to the measured outside surface temperature at midnight, or at half-hour intervals if observation data are available. This ensures the simulation starts from a physically consistent state rather than an arbitrary value.

Conduction in the ground: handled by the Force-Restore method

For the ground, the full conduction problem is more complex because the soil is a semi-infinite medium without a fixed interior temperature. Rather than solving the full 1D heat diffusion equation through the soil column, the force-restore method encodes the effect of soil conduction analytically. The restoring term $-\omega(T_{sfc} - T_{g/b})$ acts as the conduction sink, driving $T_{sfc}$ back toward the deep soil temperature $T_{g/b}$ on a timescale of $1/\omega = 86400/(2\pi) \approx 13750\ s \approx 3.8\ h$. The thermal admittance $\mu_t = \sqrt{\lambda_s \rho_s C_s}$ (where $\lambda_s$, $\rho_s$, $C_s$ are the soil thermal conductivity, density, and specific heat) controls how efficiently the soil absorbs and releases energy.

Material properties: where physics meets urban reality

The thermal conductivity $\lambda$ and thickness $e$ used in the simulation are not free parameters: they are fixed from the literature and from knowledge of the building stock. The paper assigns surface conductivity and thickness values to the four surface types present in the Toulouse center (road, wall, roof, and ground). For example, the brick walls of the Toulouse city center have higher thermal mass than glass curtain walls, which means they absorb more energy during the day and release it more slowly at night; directly contributing to the nocturnal heat island.

The sensitivity of the results to these parameters is significant: an error in $\lambda/e$ translates directly into an error in $Q_{cond}$, which then feeds back into $T_{sfc}$ through the energy balance. This is why the paper's careful material classification (including the four-category albedo scheme for wall colors) is not a cosmetic detail but a physically necessary input.

Convection transfer $Q_H$

The convective heat flux at each surface patch is parameterized through:

$$Q_H= h_f(T_a−T_{sfc})$$

The central question is: how to compute $h_f$ ? Many previous models (including SOLENE, referenced in the paper) use a constant $h_f$: a single value is applied uniformly to all surfaces. Qu et al. show this is a critical error. Instead, the model computes $h_f$ locally at every surface patch, based on the local friction velocity $u$ computed from the 3D flow solution:

$$h_f = \frac{\rho C_p u_* \kappa f_h}{σ_t \ln(\frac{d+z_0}{z_{0_T}})\sqrt{f_m}}$$

where:

• $C_p$ is the specific heat of air $(J.kg^{-1}.K^{-1}$).
• $u^∗$ is the local friction velocity $(m.s^{-1})$,computed from the $k–\varepsilon$ model.
• $\kappa$ is the von Kármán constant (0.41).
• $\sigma_t$ is the turbulent Prandtl number.
• $d (m)$ is the distance from the cell center to the wall.
• $z_0$ and $z_{0_T}$ are the momentum and thermal roughness lengths $(m)$
• $f_m$ and $f_h$ are stability functions, which reduce to 1 under neutral conditions.

This is implemented using a modified logarithmic law that accounts for atmospheric stratification through the stability functions: an explicit approach that avoids the iterative Monin-Obukhov similarity method, making it computationally tractable in 3D. The physical consequence of this choice is profound: a sun-exposed rooftop in a windy location will lose heat much more efficiently than a sheltered wall in a recirculation zone. Only a 3D flow solver can capture this variability.

The turbulence model: k–ε in the urban canopy

The authors use RANS with a k–ε turbulence closure. The choice is deliberate and honest: the k–ε model cannot resolve the large geometry-dependent eddies that LES would capture, and it overestimates dissipated energy. However, it provides acceptable accuracy at reasonable computational cost for this type of district-scale diurnal simulation. The turbulent viscosity is:

$$\mu_t=ρC_\mu \frac{k^2}{\varepsilon}$$

​ with $C_\mu=0.09$. The critical feature for urban thermal simulations is the buoyancy production term in the turbulent kinetic energy equation:

$$G_b = − \frac{\mu_t}{\rho \sigma_t} \beta g_i \frac{\partial \rho}{\partial x_i}$$

This term is what couples the temperature field to the turbulence field. When surfaces are heated, $G_b>0$, and turbulence is enhanced by buoyancy. This is particularly true near rooftop level where hot air rises from heated surfaces. This is not a correction term: in low-wind, high-sun conditions, it can be the dominant production mechanism for turbulence in the canopy.

Sensible heat flux: the turbulent engine of the urban atmosphere

The sensible heat flux at the mast level is estimated from the $k–\varepsilon$ solution as:

$$Q_H= \rho C_p w' \theta' = \rho C_p K_t \frac{\partial \theta}{\partial z}$$

where the eddy diffusivity $K_t=\frac{\nu_t}{Pr_t}$ is directly available from the turbulence model. The model-observation agreement is good during the day (differences < 40 $W.m^{-2}$) but shows a night time bias of 45–70 $W.m^{-2}$. The explanation: at night, with calm winds and a narrow surface-air temperature difference, the sensible flux is nearly zero in the observations. However, the model's initialized roof temperature is set to a measured value that creates a non-zero temperature gradient, forcing a non-zero (and here spuriously negative) heat flux.

Friction velocity: connecting turbulence to momentum

The friction velocity at the surface is defined by:

$$u^*=\left( \frac{|\tau_w|}{\rho} \right)^{1/2}$$

and at the mast level is estimated from the kinematic momentum fluxes in the two horizontal directions:

$$u^*=\left( u' w'^2+ v' w'^2 \right)^{1/4}$$

where:

$u' w' = \nu_t \frac{\partial u}{\partial z}$​ and $v'w' = \nu_t \frac{\partial v}{\partial z}$. The model captures the magnitude of $u^*$ correctly in most conditions, with difficulties appearing at the extremes (very low or very high values; a known limitation of RANS closures in complex urban geometries).

The computational domain: downtown Toulouse, not just a toy geometry

The study simulates a real urban district of 891×963×200 m centered on the intersection of Alsace-Lorraine and Pomme streets in central Toulouse, a dense 19th-century European city center with:

• Buildings 4 to 5 stories high (approximately 20 m)
• Walls primarily of brick and stone
• Roofs of clay tile
• Very sparse vegetation

The 3D urban database was provided by the city of Toulouse in AutoCAD format and preprocessed in ICEM CFD. This preprocessing was non-trivial: the raw database contained internal walls, varying heights, no ground surface, and an excess of fine detail. The preparation steps included:

• Removing internal surfaces.
• Simplifying architectural details.
• Creating the ground plane and projecting buildings onto it.
• Applying a progressive geometric detail strategy: full detail at the Alsace-Lorraine/Pomme center, simplified urban blocks in surrounding areas, and a high roughness value at the domain boundaries.

The resulting mesh is an unstructured tetrahedral grid of approximately 4.5 million cells, with resolution ranging from 0.5m near the center to 10m at the domain edges.

Boundary conditions were prescribed from measurements on the CAPITOUL meteorological mast (positioned at 47.5m above road level on a rooftop):

• Wind velocity profiles updated every 2 hours
• Wind direction updated every 2 hours
• Potential temperature profiles used to initialize the atmospheric stratification

The simulated day is 15 July 2004: a clear summer day with strong solar forcing, ideal for stress-testing the thermal model.

Results

The CAPITOUL experimental setup: what was actually measured

Before presenting the simulation results, it is worth describing what the CAPITOUL campaign actually measured, since the quality of any CFD validation is only as good as the observations it is compared against.

The meteorological mast was installed on a rooftop at 20 m height, with its top at 47.5 m above road level, at the intersection of Alsace-Lorraine and Pomme streets. It recorded continuously from February 2004 to March 2005, at one-second intervals averaged to one minute, providing:

• Short- and long-wave radiative fluxes $S^\uparrow$, $S^\downarrow$, $L^\uparrow$, $L^\downarrow$
• Sensible and latent heat fluxes $Q_H$ (via eddy-covariance)
• Wind speed, wind direction, and air temperature profiles

In addition, ten fixed infrared thermometers (IRT) were affixed to balconies and booms across three street canyons (Alsace-Lorraine, Pomme, and Rémusat), recording surface brightness temperatures of walls, roads, and roofs at fixed positions. Finally, three airborne IRT flights were performed on 15 July 2004 on board a Piper Aztec PA23 aircraft flying at 460 m altitude, producing thermal images of the district at 1.5 to 3 m resolution:

• Flight 430: 0749–0816 UTC (morning)
• Flight 431: 1115–1150 UTC (late morning)
• Flight 432: 1348–1423 UTC (early afternoon)

This dataset is exceptionally rich: it combines spatial (airborne imagery) and temporal (mast time series, fixed IRT sensors) information, which allows validation of both the spatial distribution of temperatures and their diurnal evolution.

IRT comparison: full coupling wins

The most direct test of the model is to compare the simulated brightness temperature field against the airborne IRT images. Three modelling configurations are compared against flight 431 (1138 UTC):

Model configuration	Roof temperature error	Wall error	Street error
Radiation only (no convection)	> +25°C (off scale)	Large	Large
Radiation + constant $h_f = 12\ W.m^{-2}.K^{-1}$	Overcooling ~15°C	4–5°C	4–7°C
Full radiative-convective coupling	< 10°C everywhere	< 5°C	< 3°C average

The result of the first configuration (radiation with no convective cooling) is catastrophically wrong: without $Q_H$, all absorbed radiative energy accumulates at the surface with no escape route, driving surface temperatures more than 25°C above observations. This alone justifies the coupled approach.

The second configuration (constant $h_f = 12\ W.m^{-2}.K^{-1}$) is a significant improvement, but introduces systematic overcooling at roof level of approximately 15°C. The reason is straightforward: the chosen constant value overestimates the actual spatially-averaged heat transfer coefficient at that time of day, removing too much energy from the surface. On walls (shaded or sunlit) the error is 4–5°C, and on streets 4–7°C.

Only the full radiative-convective coupling, with a spatially variable $h_f$ derived at every surface patch from the local 3D flow solution, reduces errors to below 10°C on roofs, below 5°C on walls, and below 3°C on average for streets. The remaining errors on streets are partially explained by the limited mesh resolution in the narrowest passages (a minimum of three cells per street width), which approximates shadow boundaries rather than resolving them exactly.

A key result, illustrated by the comparison between the constant and variable $h_f$ fields, is that the heat transfer coefficient varies by more than a factor of three across the domain at a given instant, depending on local wind exposure. A single constant value cannot represent this variability, regardless of how carefully it is chosen.

Diurnal evolution of brightness surface temperature $T_{br}$

The fixed IRT sensors allow validation of the temporal evolution of brightness temperature at specific positions throughout the full day of 15 July 2004. Overall, the model correctly reproduces the diurnal cycle at all sensor positions, but three instructive discrepancies are worth examining.

Alsace-Lorraine canyon walls (east and west faces): the model shows an overcooling of approximately 5°C during the evening (1800–2400 UTC) on both faces. On the west face, the model predicts a higher $T_{br}$ than observed from 0600 to 1200 UTC. This overestimation in the morning was also found by Hénon et al. (2012) using the SOLENE model on the same dataset, and is attributed to a possible underestimation by the IRT sensor itself.

Alsace-Lorraine road: the bias on the road brightness temperature can be explained by the approximation of shadow boundaries in the mesh. The value extracted from a single selected cell may differ significantly from its neighbors when the shadow edge falls nearby, making point-to-point comparison sensitive to sub-cell shadow positioning.

Pomme northeast wall: the primary disagreement occurs during the afternoon, where the bias reaches 8°C. From photographs of this wall (Figure 11 of the paper), the face is covered with numerous windows fitted with white blinds. The IRT sensor may have detected a position where blinds were closed and therefore overheated, or conversely where windows were open and interior ventilation lowered the external surface temperature. Neither effect is represented in the model, which treats the wall as a uniform opaque surface. This is a reminder that sub-facade heterogeneities in real buildings introduce irreducible noise in point-to-point validation.

Rooftop: the model shows excellent agreement before sunrise and from late afternoon to midnight. However, it exhibits a systematic advance in morning warming: the simulated roof begins heating earlier than observed. The sensor was mounted on a rooftop that, in reality, is surrounded by small obstacles (chimneys, ventilation stacks, parapets) visible in aerial photographs but absent from the simplified 3D model. These structures shade the sensor during the first hours of sunlit exposure. Once the sun has risen sufficiently to clear them, the bias disappears (consistent with the hypothesis that the error is purely geometric and not a model physics issue).

Outgoing radiative fluxes: $L^\uparrow$ and $S^\uparrow$

The meteorological mast measures both the outgoing long-wave flux $L^\uparrow$ and the outgoing short-wave (solar) flux $S^\uparrow$, providing a direct check on the radiative model.

Long-wave flux $L^\uparrow$: the outgoing infrared flux is directly tied to surface temperature through $L^\uparrow = \varepsilon \sigma T_{sfc}^4$. Because of the fourth-power dependence, errors in $T_{sfc}$ are amplified in $L^\uparrow$: a 3°C error at 320 K corresponds to approximately $20\ W.m^{-2}$. The model shows:

• An underestimation of approximately $20\ W.m^{-2}$ during the night, directly traceable to the nighttime surface temperature bias discussed above.
• An overestimation of up to $100\ W.m^{-2}$ at noon, traceable to the morning advance in roof warming.

The difference between the modeled $L^\uparrow$ computed at the roof surface and at the mast level is generally less than $20\ W.m^{-2}$, confirming that the vertical variation of the outgoing longwave flux between the surface and 47.5 m is small in a transparent atmosphere.

Beyond the scalar comparison at the mast, the model allows full 3D visualization of the $L^\uparrow$ distribution across the domain (Figure 14 of the paper, first figure shown in this article). This reveals significant spatial variability: horizontal surfaces (rooftops, roads) are substantially warmer than vertical surfaces (walls) in the daytime, and therefore emit significantly more longwave radiation. This cannot be captured by any model that represents the urban surface as a single averaged layer.

Short-wave flux $S^\uparrow$: For the gray roof at the mast position, an albedo of $\alpha = 0.15$ was assigned, consistent with the literature. The agreement between model and observation is very good throughout the day, with an overestimation of only approximately $20\ W.m^{-2}$ at noon. This validates both the solar radiation scheme and the albedo classification strategy. The 3D visualization of $S^\uparrow$ (Figure 16 of the paper) also makes immediately visible which surfaces have higher albedo (brighter in the image) and the sharp shadow boundaries cast by building edges — both direct outputs of the DOM solver on the full 3D mesh.

Sensible heat flux $Q_H$ and friction velocity $u^*$

Sensible heat flux: the sensible heat flux at the roof surface is compared to the eddy-covariance measurement on the mast. The agreement is good during the day, with differences generally below $40\ W.m^{-2}$. The main discrepancy occurs at night: from midnight to 0700 UTC, the observed $Q_H$ is nearly zero (the air just above the roof is nearly at the same temperature as the roof surface itself). However, the model initializes the roof temperature to a measured value that is lower than the air temperature at that hour, producing a spuriously negative sensible heat flux (i.e., heat flowing from the air toward the surface). The nighttime bias of 45–70 $W.m^{-2}$ is therefore a consequence of the initialization procedure rather than a model physics failure. From 1300 to 1500 UTC, a secondary discrepancy arises that the authors attribute to the influence of unresolved rooftop geometry on local convective exchange.

Friction velocity: the friction velocity $u^*$ provides a check on the momentum transfer, independent of the thermal model. The model captures the correct order of magnitude throughout the day, with good agreement when $u^* < 0.2\ m.s^{-1}$. At higher values, the model tends to overestimate $u^*$ at the mast level during the night (1800–2400 UTC). This is a known limitation of the RANS $k$–$\varepsilon$ closure in complex geometries: it overestimates turbulent momentum transfer in separated and recirculating flows, and estimates of $u^*$ from sonic anemometry are themselves uncertain when values are small. The difference between $u^*$ computed at the roof surface and at the mast level is also partially a numerical artifact of the gradient estimation method.

The critical application: how thermal effects reshape pollutant dispersion

Having established confidence in the thermal model through the CAPITOUL comparisons, the authors use the same framework to study a question of direct public health relevance: what happens to a pollutant released at ground level in a city center on a hot summer day, compared to the standard neutral-atmosphere assumption?

Two simulations are performed with identical wind conditions ($U_{ref} = 3\ m.s^{-1}$ westward at 47.5 m, neutral temperature profile at 22°C), differing only in the treatment of thermal effects:

• Neutral case: no thermal effects, no surface heating, no buoyancy.
• Thermo-radiative case: full 3D radiative-convective-conductive coupling active, solar forcing from 1000 UTC for 3 hours.

A passive tracer is released at ground level at the center of the domain for 3 hours.

Effect on turbulent kinetic energy: in the neutral case, TKE production is purely mechanical, confined to the immediate vicinity of building edges where flow separates. With thermal coupling active, heated surfaces continuously inject buoyant energy into the flow. TKE increases throughout the entire domain, not just near buildings, and is particularly enhanced near rooftop level where the buoyancy production term $G_b$ is largest. The flow never reaches a stationary state in the thermo-radiative case: the continuously evolving surface temperatures continuously modify the buoyancy forcing, making the flow inherently unsteady even under constant wind boundary conditions.

Effect on ground-level concentration: in the neutral case, the plume travels westward with the dominant wind. Inside courtyards and narrow street canyons, wind speeds are dramatically reduced, creating stagnation zones where the pollutant accumulates and is poorly diluted. High ground-level concentrations persist well downwind of the source.

With thermal transfer active, surface heating drives a persistent upward buoyant motion. The plume rises from the point of emission, concentrations near the ground are lower downstream, and the courtyard concentrations are reduced. However, concentrations above roof level are higher: the thermal uplift that protects pedestrians at street level transfers the pollution upward into the layer where it can be advected over longer distances.

These are not subtle quantitative differences. They represent qualitatively different dispersion regimes: the neutral model predicts a ground-hugging plume trapped in street canyons, while the thermally coupled model predicts a rising plume escaping the canopy. If you are a city air quality manager assessing the risk of a ground-level accidental release on a hot summer afternoon, the neutral assumption gives you the wrong geometry of the hazard zone: not just wrong concentrations, but wrong spatial distribution, wrong affected population, and wrong mitigation strategy.

From simulation to design: what this enables for urban planners

The CAPITOUL validation gives confidence that Code_Saturne's thermo-radiative-convective framework can be trusted as a digital laboratory for urban design interventions. The authors themselves enumerate the applications:

🏗️ High-Albedo Surfaces ("White Roofs")

In the wall thermal model, changing the albedo α directly modifies the absorbed shortwave radiation:

$$(1−\alpha)(S_D+S_f+S_e)$$

Increasing α from 0.15 (gray roof, as measured in Toulouse) to 0.7 (white roof coating) reduces the solar energy absorbed by the surface by 65%. This has a cascade effect: lower surface temperature → lower $Q_H$ into the street air (daytime cooling) and lower $Q_cond$ stored in the wall (reduced nocturnal heat release) → cooler air temperature both day and night. The model can quantify this cascade for a specific district with its specific geometry and orientation.

🌿 Green Roofs and Green Areas

The framework can be extended to include the latent heat flux term (currently neglected) to represent vegetation. Evapotranspiration removes energy from the surface as latent heat rather than sensible heat, directly cooling the local air. The model provides the spatial wind and turbulence fields needed to correctly compute where this cooling has the most impact which is a function of local flow patterns that vary enormously within a single city block.

🌬️ Street Canyon Ventilation

By running the full 3D flow simulation, the model identifies which streets are well-ventilated (high local $u^∗$, high $h_f$, efficient heat removal) and which are thermal traps (low wind speed, low $h_f$, accumulating heat). This information can directly inform decisions about building height ratios, street widths, and building orientation in new developments.

☁️ Air Quality "Hot Spots"

As the pollutant dispersion results show, thermal effects fundamentally change where pollution concentrates. The coupled model can identify specific locations in a city center where ground-level emissions are most likely to accumulate under typical summer heatwave conditions, and evaluate whether a cool pavement or a row of trees would meaningfully change the local air quality.

Why code_saturne? A reference for open-source scientific CFD

code_saturne is not a commercial black box. Developed and maintained by EDF R&D, it is:

✅ Open-source (GNU GPL license): every equation described in this article is inspectable, modifiable, and fully auditable

✅ Finite-volume on unstructured meshes: handles the extremely complex geometries of real city centers without simplification

✅ Physically complete: the atmospheric module includes RANS (k–$\varepsilon$) and LES turbulence, DOM radiation (short-wave and long-wave separately), the force-restore ground model, the wall thermal model, and passive scalar dispersion, all fully coupled

✅ Massively parallel: scales to thousands of CPU cores, making district-scale diurnal simulations tractable

✅ Validated for real urban environments: the CAPITOUL study is the most comprehensive validation to date, but it builds on prior validation against the MUST field campaign and idealized cases

✅ Free to use: which means cities of all sizes, developing countries, and academic research groups can access the same tools used by EDF and CEREA

Getting started: your first urban thermal simulation

Here is a practical roadmap to replicate this type of study for your own city:

• Geometry: Obtain the 3D urban database from your municipality (many European cities have this in CityGML or similar format). Import it into ICEM CFD or Salome-Meca (open-source). Perform geometry cleanup: remove internal walls, create the ground surface, project buildings, and simplify detail progressively from domain edge to center of interest.

• Mesh: Build an unstructured tetrahedral mesh. Target ~0.5 m resolution at your area of interest, coarsening to ~10 m at the domain boundaries. Total cell count of 3–5 million is a practical starting point for a district of ∼1 km².

• Material Properties: Classify your building surfaces by albedo and emissivity based on visual inspection or municipal records. At minimum, distinguish roof materials (tile, gravel, bitumen) and wall colors. Define wall thermal conductivity λ and thickness e from the literature for your building type: these control the conductive thermal resistance $R=e/\lambda$ and therefore the building's thermal inertia.

• Atmospheric Boundary Conditions: Use measured wind speed, direction, and temperature profiles from the nearest meteorological station, updated every 2 hours. Set the atmospheric stability from potential temperature profiles.

• Activate Full Coupling: In Code_Saturne's atmospheric module, activate the DOM radiative transfer solver (32 directions for a first run, 128 for validation), the force-restore ground model, and the wall thermal model. Set the solar position from your geographic coordinates and the simulation date.

• Validate: Compare with any available surface temperature measurements, ideally including a rooftop meteorological station. If possible, use satellite or drone-based thermal infrared imagery for spatial validation. Remember to compare brightness temperatures, not raw surface temperatures.

• Design Exploration: Once validated, change one variable at a time (roof albedo, surface material, tree canopy, building height) and quantify the effect on surface temperature, sensible heat flux, and pollutant dispersion. This is your digital urban laboratory.

Conclusion: The city as a solvable equation

The work of Qu, Milliez, Musson-Genon, and Carissimo is more than an academic study. It is a proof of concept that the urban thermal environment is quantitatively predictable in three dimensions, at the building-resolving scale, over a full diurnal cycle, in a real European city center.

The physics is encoded in the surface energy balance. The radiative exchange is encoded in the Discrete Ordinates Method. The surface temperature emerges from the coupled force-restore and wall thermal models. Wall conduction links the exterior surface temperature to the building's thermal mass, encoding the slow nocturnal heat release that sustains the urban heat island long after sunset. The convective cooling is computed locally at every surface patch from the 3D flow solution. And Code_Saturne ties them all together into a solvable, reproducible, open-source simulation framework.

Perhaps the most sobering finding of the paper is the pollutant dispersion result: assuming a neutral atmosphere (as the vast majority of urban flow studies still do today) gives you the wrong plume shape, the wrong ground-level concentrations, and the wrong identification of at-risk locations. On a hot summer day in a city center, buoyancy is not a perturbation. It is a first-order effect.

We have the tools. We have the physics. We have the open-source software. We even have the validation data from a year-long field campaign in Toulouse. What we need now is more engineers, planners, and researchers willing to run the simulations and bring quantitative, physics-based evidence into the rooms where urban design decisions are made.

Because when the next heatwave kills people in city centers (and it will) the question we will have to answer is not "did we have the tools to understand this?" (we clearly do) but "did we use them?"

References

📌 The full paper: Qu Y., Milliez M., Musson-Genon L., Carissimo B. — "Modelling Radiative and Convective Thermal Exchanges over a European City Center and Their Effects on Atmospheric Dispersion" — Sustainability, 2022, 14(12), 7295. DOI: 10.3390/su14127295

📌 The CAPITOUL experimental dataset: Masson V. et al. — "The Canopy and Aerosol Particles Interactions in TOulouse Urban Layer (CAPITOUL) Experiment" — Meteor. Atmos. Phys., 2008, 102, 135–157.

[1]: Butera, Federico. (2019). ENERGY AND RESOURCE EFFICIENT URBAN NEIGHBOURHOOD DESIGN PRINCIPLES FOR TROPICAL COUNTRIES - A Practitioner's Guidebook.

[2]: Pigeon, G., Legain, D., Durand, P., & Masson, V. (2006). Anthropogenic heat release in an old European agglomeration (Toulouse, France). International journal of climatology, 27(14), 1969-1981.