A high-resolution urban CO₂ transport model with anthropogenic and biogenic fluxes

We develop PALM-CO₂, a high-resolution urban carbon dioxide transport model with anthropogenic and biogenic carbon emissions. The PALM-CO₂ model is implemented within the open-source urban-flow large-eddy simulation (LES) model PALM. Anthropogenic CO₂ emissions are prescribed using an external emission inventory, while biogenic CO₂ fluxes are computed online using the Vegetation Photosynthesis and Respiration Model (VPRM). In addition, custom output modules are developed to diagnose carbon fluxes. PALM-CO₂ is validated through a case study in London, comprising an 8 by 8 km² domain covering the borough of Camden at a resolution of 10 m. Simulations are driven by reanalysis meteorological forcing and background CO₂ concentrations, while the hourly anthropogenic emissions at 10 m resolution are explicitly derived in this study. Validation against eddy-covariance flux measurements within the study region confirms that the model captures the diurnal variation of turbulent transport and anthropogenic emissions. The biogenic flux module is evaluated independently using observed monthly diurnal profiles of biogenic CO₂ fluxes from a deciduous forest site in the Czech Republic. The comparison shows that the model captures realistic seasonal and diurnal variations in ecosystem carbon exchange. Analysis of the London simulations reveals strong spatial heterogeneity in near-surface CO₂ concentrations arising from the interaction of building-induced turbulence, diurnal boundary-layer evolution, and emission patterns. PALM-CO₂ provides a high-resolution framework for investigating CO₂ transport processes in complex urban and vegetated environments, providing improved quantification of urban emission sources.

1 Introduction

Greenhouse Gas (GHG) emissions cause global climate change, with carbon dioxide (CO₂) being the most significant contributor (Intergovernmental Panel on Climate Change, 2023). To understand the impact of CO₂ emissions, physical models have been used for modelling the carbon cycle and climate. The Coupled Model Intercomparison Project (CMIP) (Anav et al., 2013) compares 18 earth system models' ability to simulate the global carbon cycle across land and ocean; the models provide important insights for stakeholders, for example, climate projections from these models are used in Intergovernmental Panel on Climate Change (IPCC) report.

In the global carbon budget (Friedlingstein et al., 2023) for 2013–2022, anthropogenic activities (fossil CO₂ and land-use change) account for 12.9 GtC per year, while land uptake is 3.3 GtC per year and ocean uptake is 2.8 GtC per year. The numbers show the imbalance of the emission caused by anthropogenic activities and the huge uptake provided by biospheres and natural processes. Accurately representing both anthropogenic emissions and biospheric CO₂ fluxes is therefore essential for reliable carbon dioxide modelling.

A large portion of emissions occurs in urban environments (Dhakal, 2010; Marcotullio et al., 2013) due to industrialisation and high population density – the energy sector produces 40 %–80 % of global energy-related carbon emissions. Urban environments exhibit strong spatial and temporal heterogeneity in emissions as a result of diverse land use, building topography, and socioeconomic activities, leading to complex intra-urban emission patterns. For example, Christen et al. (2011) estimated carbon emissions for a Vancouver neighbourhood at a 50 m grid resolution using activity data from buildings, waste, and transportation, demonstrating how heterogeneous urban land use influences emission patterns. Smith et al. (2025) adopted land-use regression to model carbon dioxide mixing ratios in the San Francisco Bay Area, showing the great intra-urban variance of carbon emission patterns. These studies have indicated the complexity in spatial and temporal patterns of anthropogenic emissions.

For urban carbon transport modelling, anthropogenic emission inputs are typically derived from inventories at kilometre-scale or coarser resolution and must be spatially disaggregated to match the model grid. This downscaling process commonly employs surrogate variables such as population density, building data, traffic activity, and land-use. While the methodology for generating fine-scale emission maps is acknowledged, their implementation depends strongly on the availability of local activity and geospatial datasets. Several urban inventories have achieved resolutions from 1 km down to 50 m by integrating activity data and remote sensing products (London Datastore, 2023; Christen et al., 2011; Järvi et al., 2019; Cai et al., 2023). However, 10 m resolution anthropogenic CO₂ emission datasets suitable for large-eddy simulation applications remain scarce. This is a critical limitation, as coarse-resolution inventories cannot resolve the strong spatial heterogeneity of urban emissions, leading to biased or overly smoothed surface flux representations at neighbourhood scales.

In addition to anthropogenic emission, biogenic fluxes can serve as a source or sink for atmospheric CO₂ (Friedlingstein et al., 2023). The land biosphere and oceans together currently act as a net sink of anthropogenic CO₂, absorbing about half of human emissions each year (Friedlingstein et al., 2023). However, on the urban scale, the diurnal cycle of biosphere activities yields temporally varying biogenic fluxes. To represent these processes, biogenic carbon models have been developed that explicitly simulate the underlying biochemical mechanisms of vegetation photosynthesis and ecosystem respiration. For example, Gerbig and Koch (2024) applied the Vegetation Photosynthesis and Respiration Model (VPRM, Mahadevan et al., 2008) to calculate biosphere-atmosphere exchange fluxes for CO₂ across Europe. At the urban scale, Järvi et al. (2019) introduced a new biogenic carbon model that accounts for both human metabolism and vegetation effects, while Stagakis et al. (2025) conducted an intercomparison of four different biogenic models, highlighting discrepancies between model formulations. These studies demonstrate that biogenic fluxes can substantially influence urban CO₂ budgets, yet their integration with high-resolution transport models remains limited.

To investigate the transport of atmospheric components, physically based urban flow models have been developed to resolve local meteorological conditions. These models are grounded in physical principles and solved numerically, enabling simulations at very fine spatial-temporal resolutions. The models have been widely used to study heatwave (Anders and Maronga, 2025) and pollution episodes (An et al., 2007), and local micro-climate's impacts on, for example, building energy usage (Pfafferott et al., 2021) and safety operations of air vehicles (Jiang et al., 2024) as well as urban ecosystems (e.g. vegetation degradation (Neves et al., 2018)). Such models generally take consideration of local urban topology including terrain height and building geometries, land and vegetation types, meteorological driving forces (e.g. wind and radiation), as well as anthropogenic activities (e.g. chemical emissions and heat fluxes). Examples of such models are PALM (Maronga et al., 2020) and DALES (Heus et al., 2010).

Beyond advancing the understanding of physical mechanisms, finer-resolution carbon transport models in urban environments are essential for accurately representing concentration fields, thereby improving the model’s capability to infer urban emission sources from ground-based measurements. This has been widely acknowledged within air pollution studies, as fine-resolution modelling enhances exposure assessment and supports source attribution and mitigation strategies (Stroh et al., 2007; Li et al., 2016; Fenech et al., 2018). In contrast, such approaches are less commonly applied in the CO₂ study. For example, Brunner et al. (2019) demonstrated that neglecting realistic vertical emission profiles can lead to substantial overestimation of near surface CO₂ concentrations by up to 14 % in summer and 43 % in winter – highlighting the importance of vertical resolution in modelling. The limited use of fine scale CO₂ transport modelling in urban environments may partly stem from the scarcity of dense observational datasets. This gap motivates the development of integrated, high-resolution modelling frameworks that consistently couple urban flow dynamics with both high spatial resolution anthropogenic CO₂ emission datasets and biogenic CO₂ fluxes.

In this study, we develop an urban CO₂ transport modelling framework (PALM-CO₂¹) with 10 m resolution in both the horizontal and vertical directions, coupled with a CO₂ emission inventory resolved at the same spatial scale. The model is built upon the PALM large eddy simulation (LES) system (Maronga et al., 2020), which provides an accurate representation of urban scale flows and turbulence. Within this framework, we implement a biogenic CO₂ flux module and incorporate anthropogenic emission inventories as external inputs, enabling the simultaneous simulation of transport processes and surface–atmosphere CO₂ exchange in complex urban environments.

A key component of this work is the development of a new 10 m anthropogenic CO₂ emission inventory derived from an existing 1 km inventory using multiple spatial proxies. This high-resolution inventory is combined with high-resolution biogenic CO₂ fluxes and implemented as spatially and temporally varying surface scalar fluxes within the PALM LES framework. The resulting modelling system enables the representation of neighbourhood-scale variability in both anthropogenic emissions and biogenic uptake, which cannot be captured using coarse-resolution inventories alone.

To support urban-scale applications, we design a dedicated workflow to spatially disaggregate anthropogenic CO₂ emissions for the study area of London. This workflow generates a high resolution emission grid map at 10 m × 10 m resolution by integrating multiple sources of locally available urban geoinformation, including residential and commercial property datasets, gas consumption statistics, industrial point-source inventories, and detailed road network maps. These datasets are combined to apportion emissions consistently across urban surface types and activity sectors, ensuring compatibility with the PALM computational grid.

Model performance is evaluated against eddy-covariance measurements within the study area. The results show that PALM-CO₂ reproduces the observed diurnal variability of turbulent CO₂ transport and anthropogenic emission signals. The biogenic flux module is evaluated independently using observed monthly diurnal CO₂ fluxes from a deciduous forest site in the Czech Republic, demonstrating its ability to capture realistic seasonal and diurnal patterns of ecosystem carbon exchange.

The London simulations reveal pronounced spatial heterogeneity in near-surface CO₂ concentrations arising from the combined effects of urban morphology, building-induced turbulence, boundary-layer evolution, and spatially varying emission sources. PALM-CO₂ provides a high-resolution framework for investigating urban CO₂ transport and surface-atmosphere exchange processes and offers new opportunities for urban carbon monitoring, source attribution, and future emission inversion studies.

2 Model description: PALM-CO₂

PALM-CO₂ is implemented within the PALM model system (Maronga et al., 2020), a large-eddy simulation (LES) framework designed for high-resolution simulations of atmospheric flow in complex urban environments. PALM explicitly resolves the dominant turbulent eddies while parameterizing subgrid-scale (SGS) motions, and provides a modular structure that supports detailed representations of urban topography, land–atmosphere interactions, and scalar transport. Building upon these capabilities, PALM-CO₂ extends the existing trace-gas functionality of PALM to enable the simulation of urban-scale CO₂ transport coupled with spatially and temporally resolved surface fluxes.

2.1 PALM framework and LES formulation

In PALM, the transport of passive scalar quantities, including trace gases, is governed by the filtered scalar conservation equation (Maronga et al., 2020):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{s}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \mathit{\rho }{\stackrel{\mathrm{‾}}{u}}_j\stackrel{\mathrm{‾}}{s}}{\partial x_j}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial }{\partial x_j}}\left(\mathit{\rho }{\stackrel{\mathrm{‾}}{u}}_j^{\prime \prime }{s}^{\prime \prime }\right)+{\mathit{\chi }}_{\mathrm{s}},\end{array}$$

where s is the concentration of the trace gas in kg m⁻³, t is time, ρ is the density of dry air, u_j is the jth component of velocity, x_j is the spatial coordinate, and χ_s is the source or sink term of the passive scalar. Variables with double prime (^′′) refer to the subgrid-scale (SGS) component of the corresponding variables, while the over-bar indicates the filtered quantities. The advection and SGS turbulent flux divergence terms are discretized using finite-difference schemes implemented in PALM, and SGS fluxes are parametrised via an eddy-diffusivity approach consistent with the LES formulation.

PALM-CO₂ leverages these existing scalar transport codes and augments them with dedicated modules to represent anthropogenic and biogenic CO₂ fluxes. These fluxes enter the scalar transport equation exclusively through the source or sink term χ_s, ensuring a clean separation between atmospheric transport processes and surface exchange parametrisations.

Figure 1 illustrates the overall structure of PALM-CO₂, including the coupling between the flow solver, the biogenic flux module, and the prescribed anthropogenic emission inputs.

Figure 1 Schematic workflow of the PALM–VPRM modelling system. Blue boxes indicate static and dynamic drivers, green boxes indicate VPRM-related vegetation inputs, purple boxes denote model components, orange boxes represent anthropogenic emissions inputs, and black boxes indicate model outputs and post-processing. Arrows show the direction of data exchange and processing. Input files are indicated in a monospace font.

2.2 Biogenic flux model

Biogenic CO₂ exchange is represented using the Vegetation Photosynthesis and Respiration Model (VPRM; Mahadevan et al., 2008), which parametrises the ecosystem-scale carbon uptake and release based on vegetation state and meteorological drivers. VPRM has been widely applied in regional and urban-scale carbon studies (Ahmadov et al., 2007, 2009; Pillai et al., 2010, 2011). The parameters in the VPRM model are typically fitted using measurements from local eddy-covariance flux towers. In this study, we select a set of parameters (see Appendix F) from the WRF-Chem module (Skamarock et al., 2019) that has been calibrated using flux tower measurement in Europe (Lian et al., 2021).

In PALM-CO₂, VPRM computes gross primary production (GPP) and ecosystem respiration (RESP) using the following equations:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\text{GPP}& =\mathit{\lambda }\times T_{\text{scale}}\times P_{\text{scale}}\times W_{\text{scale}}\\ & \times \text{EVI}\times {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\text{PAR}/{\text{PAR}}_{\mathrm{0}}}}\times \text{PAR},\end{array}\end{array}$$

where λ is a parameter fitted with the measurement of biogenic fluxes. Vegetation status and water stress are represented by the Enhanced Vegetation Index (EVI) and the Land Surface Water Index (LSWI). P_scale is a coefficient accounting for the degree of leafiness which is computed from LSWI and also depends on EVI to determine whether the vegetation has reached full leaf expansion,

$$\begin{array}{}\text{(3)}& P_{\text{scale}}=\left\{\begin{array}{ll}\mathrm{1},& \left(\text{EVI}>{\text{TH}}_{\text{leaf full expansion}}\right),\\ \frac{\mathrm{1}+\text{LSWI}}{\mathrm{2}},& \left(\text{EVI}\le {\text{TH}}_{\text{leaf full expansion}}\right).\end{array}\right.\end{array}$$

TH_{leaf full expansion} is a threshold for vegetation growing cycle. EVI above the threshold means the growing period between leaf full expansion and senescence and below the threshold is during bud burst to leaf full expansion and during senescence. The threshold is computed from the annual maximum and minimum EVI for the computation grid,

$$\begin{array}{}\text{(4)}& {\text{TH}}_{\text{leaf full expansion}}={\text{EVI}}_{\text{min}}+\mathrm{0.55}\times ({\text{EVI}}_{\text{max}}-{\text{EVI}}_{\text{min}}).\end{array}$$

W_scale accounts for the water content in the vegetation and is computed from LSWI,

$$\begin{array}{}\text{(5)}& W_{\text{scale}}=\left\{\begin{array}{ll}\frac{\text{LSWI}-{\text{LSWI}}_{\text{min}}}{{\text{LSWI}}_{\text{max}}-{\text{LSWI}}_{\text{min}}},& \left(\text{grassland and savanna,}\right),\\ \frac{\mathrm{1}+\text{LSWI}}{\mathrm{1}+{\text{LSWI}}_{\text{max}}},& \left(\text{all other classes}\right).\end{array}\right.\end{array}$$

PAR is photosynthetically active radiation which is correlated with short-wave (SW) radiation: $\text{PAR}\approx \text{SW}/\mathrm{0.505}$; PAR₀ is a parameter for half photosynthetic saturation. T_scale is a coefficient computed from air temperature T,

$$\begin{array}{}\text{(6)}& T_{\text{scale}}={\displaystyle \frac{(T_{\text{air}}-T_{\text{min}})(T_{\text{air}}-T_{\text{max}})}{(T_{\text{air}}-T_{\text{min}})(T_{\text{air}}-T_{\text{max}})-(T_{\text{air}}-T_{\text{opt}}{)}^{\mathrm{2}}}},\end{array}$$

with T_min, T_opt, and T_max denoting minimum, optimal, and maximum temperatures for photosynthesis. T_air is the air temperature.

The respiration flux is mainly determined by the air temperature T_air and is computed from

$$\begin{array}{}\text{(7)}& \text{RESP}=\mathit{\alpha }{T}_{\text{air}}+\mathit{\beta },\end{array}$$

where α and β are coefficients that vary among different vegetation types. The net ecosystem exchange (biogenic flux) is then calculated

$$\begin{array}{}\text{(8)}& \text{NEE}=-\text{GPP}+\text{RESP},\end{array}$$

where a negative value means net uptake and a positive value means net emission.

2.3 Diagnostic turbulent CO₂ fluxes

To facilitate comparison with eddy-covariance observations and to analyse vertical CO₂ exchange, PALM-CO₂ includes a diagnostic computation of turbulent scalar fluxes. The vertical turbulent flux (F) is defined as

$$\begin{array}{}\text{(9)}& F=\langle w^{\prime }{s}^{\prime }\rangle ,\end{array}$$

where angle brackets denote the temporal average,

$$\begin{array}{}\text{(10)}& \langle \mathit{\varphi }\rangle ={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+\mathit{\tau }}{\int }}\mathit{\varphi }\left(t\right)\mathrm{d}t,\end{array}$$

where τ is the time window for averaging. w^′ and s^′ are the deviations from the mean of the vertical velocity component and tracer gas concentration (e.g. CO₂) respectively,

$$\begin{array}{}\text{(11a)}& {\displaystyle }{w}^{\prime }=w-\langle w\rangle ,\text{(11b)}& {\displaystyle }{s}^{\prime }=s-\langle s\rangle .\end{array}$$

Then the flux can be computed from the following temporal averaged quantities,

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\displaystyle }\langle w^{\prime }{s}^{\prime }\rangle & {\displaystyle }=\langle (w-\langle w\rangle )(s-\langle s\rangle )\rangle \\ {\displaystyle }& {\displaystyle }=\langle ws-\langle w\rangle s-w\langle s\rangle +\langle w\rangle \langle s\rangle \rangle \\ {\displaystyle }& {\displaystyle }=\langle ws\rangle -\langle w\rangle \langle s\rangle \\ {\displaystyle }& {\displaystyle }=\langle \stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{s}\rangle -\langle \stackrel{\mathrm{‾}}{w}\rangle \langle \stackrel{\mathrm{‾}}{s}\rangle .\end{array}\end{array}$$

In the last line, the spatially filtered quantities $\stackrel{\mathrm{‾}}{w}$ and $\stackrel{\mathrm{‾}}{s}$ and their product $\stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{s}$ from the Large-eddy simulation (LES) code are used to compute the vertical turbulent flux for the CO₂ concentration. A major assumption is the omission of the sub-grid scale turbulent flux, which requires the grid size to be small enough to resolve the total turbulent flux. It is shown in the grid sensitivity study (Sect. 4.4) that the resolved flux at medium height boundary layer (higher than three roughness-layer heights) takes over 96 % of the total turbulent flux. Thus, the assumption is valid. A minor assumption is that the spatial filtering and time averaging commute and that the filter preserves the mean, such as $\langle \cdot \rangle =\langle \stackrel{\mathrm{‾}}{\cdot }\rangle$. This assumption requires the averaging window is long enough to converge the statistics at all spatial locations. A computation of eddy turnover time (about 450 s for the grid sensitivity study case) shows that a half-hour average window (1800 s) covers about four eddy turnover periods, validating this minor assumption. The temporal averaging of $\stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{s}$ is added to the code via a user-defined procedure. Its rolling sum is computed during time-stepping in the code, and its temporal average 〈ws〉 is computed and reported after each averaging time window τ. 〈w〉 and 〈s〉 can be easily recorded using the default output settings in PALM. The time window for averaging is set as 0.5 h.

2.4 Anthropogenic emission representation and disaggregation methodology

Anthropogenic CO₂ emissions for urban areas are commonly provided as coarse-resolution inventories, typically with a kilometre-scale spatial resolution, and as annual totals. For use in high-resolution urban transport models, such inventories must be disaggregated both spatially and temporally to match the computational grid and time step. A direct interpolation of coarse grid values onto a fine-resolution model grid is straightforward but insufficient to represent realistic urban emission patterns, which are often dominated by localised sources such as building heating systems, traffic corridors, and industrial point sources. The method for processing emission inputs for chemistry transport models has been implemented in several existing tools, including the online emission module in COSMO (Jähn et al., 2020) and the stand-alone processors HERMES (Guevara et al., 2019), EMIPS (Chen et al., 2023), and FUME (Belda et al., 2024). However, these tools are primarily designed for global and regional applications and directly use emission inventories at those scales; therefore, they are not readily applicable to our study.

To address this, a sector-based disaggregation strategy is adopted. Disaggregation is applied on the sub-subsection levels defined in the LAEI dataset (London Datastore, 2023, listed in Appendix A). Emission sub-subsectors contributing only a small fraction of total emissions are distributed uniformly in space and time after interpolation to the fine grid. In contrast, dominant sub-subsectors are spatially disaggregated using additional geospatial surrogate data that represent the underlying emission-generating activities. This approach balances accuracy with data availability, as comprehensive surrogate datasets are not uniformly available for all emission sub-subsectors. The selection of sub-subsectors for detailed disaggregation is therefore based on their relative contribution to total emissions.

For sub-subsectors selected for fine-scale treatment, spatial disaggregation relies on auxiliary geospatial datasets such as building-level energy information, industrial point-source inventories, and road network data. These datasets are used as spatial proxies to redistribute emissions from coarse administrative or grid-based units to individual buildings, road segments, or point sources, which are then aggregated to the target model grid resolution.

Temporal disaggregation is applied to convert annual emissions into hourly fluxes. Sub-subsector-specific temporal profiles are used where available, capturing diurnal, weekly, and seasonal variability in anthropogenic activities. Sub-subsectors with negligible contributions are assigned constant temporal profiles, while dominant sub-subsectors employ time-varying profiles derived from established emission databases.

2.5 Implementation of PALM-CO₂ framework

PALM provides an interface for user-defined chemical mechanisms where one can set up the chemical components (Khan et al., 2021). In PALM-CO₂, CO₂ is represented as a single passive scalar without chemical reactions. Both anthropogenic and biogenic fluxes are treated as surface source terms and are incorporated into the scalar transport equation through χ_s in Eq. (1).

Biogenic CO₂ fluxes are computed by VPRM, whereas anthropogenic emissions are provided as externally prepared, time-resolved gridded datasets. VPRM (described in Sect. 2.2) is implemented in PALM as part of an online chemistry module, where biogenic fluxes are computed prior to solving the scalar transport equation using near-surface air temperature and short-wave radiation retrieved from PALM runtime variables. Anthropogenic emission sources (described in Sect. 3.2) are supplied as hourly NetCDF files, requiring activation of the level-of-detail (LOD) 2 option in PALM.

Both biogenic and anthropogenic fluxes are prescribed as source terms in the lowest atmospheric grid cell adjacent to the surface. In locations where complex terrain or building structures elevate the first fluid grid cell above ground level, the surface fluxes are applied at that elevated cell. This approach maintains a physically consistent coupling between surface emissions and the resolved airflow, ensuring that fluxes are introduced directly into the urban canopy layer where turbulent exchange and transport are explicitly represented.

3 Case Study and Experimental Setup

3.1 Study area

A study area in London was selected for model validation and analysis, as shown in Fig. 2. London is selected for its high population density, high level of urbanisation, and its vast amount of green infrastructures. In Camden, there is a CO₂ flux measurement station located at 190 m height on the BT tower (Helfter et al., 2011, 2016), providing an independent dataset for evaluating the performance of the model developed in this study.

Figure 2 A nested study area surrounding M25 motorway and London Borough of Camden. Coordinates are in British National Grid coordinate system. (a) The parent domain (denoted as N01) is roughly surrounding M25 motorway and the Greater London Area (GLA); the domain size is 60 km by 60 km. (b) The child domain (denoted as N02) encompasses Camden and its size is 8 km by 8 km. A measurement site (BT tower) is marked with a green star in the child domain. A x=4000 m line where the vertical cross-section data is extracted for analysis is shown in blue; the origin of the N02 domain is the bottom-left corner. Map data in the figure are from OpenStreetMap (© OpenStreetMap Distributed under the Open Data Commons Open Database License (ODbL) v1.0 (for further information, please see OpenStreetMap)).

The Greater London Authority (GLA) area (Fig. 2a) is characterised by a highly heterogeneous land-use and urban landscape. At its core lies a dense, multifunctional urban structure dominated by commercial, financial, and high-density residential uses, with extensive transport infrastructure and limited open space. Moving outward, land use transitions to predominantly residential neighbourhoods of medium to low density, interspersed with local centres, industrial estates, and major road and rail corridors. Despite its metropolitan character, the GLA area contains a substantial proportion of green and blue spaces, including large public parks, river corridors – most notably the River Thames – and protected areas (Wilby and Perry, 2006).

Sitting at the core of the GLA, the borough of Camden (Fig. 2b) reflects this heterogeneous landscape. From north to south, there is a mixing of green spaces (Hampstead Heath and the Regent's Park), residential area, and commercial and business centres around three main train stations (Euston, St Pancras and King's Cross). These are also shown in Fig. 3 which depicts a bird's-eye view of dense building areas with parks and stations. Overall, land use in the GLA and Camden shows a complex mosaic in which intense urban development coexists with natural and semi-natural landscapes, resulting in strong spatial contrasts across relatively short distances.

A nested domain is used, which consists of a coarse-resolution domain (parent domain, denoted as N01 domain) and a fine-resolution domain (child domain, denoted as N02 domain). The parent domain (GLA)'s range is 60 km by 60 km in the horizontal and 3 km in the vertical, while the grid resolution in $x,y,z$ is 100, 100, and 50 m respectively. The child domain (Camden)'s range is 8 km by 8 km in the horizontal and 1 km in the vertical, while the grid resolution is 10 m in all directions.

Figure 3 Bird's-eye view of buildings located in a region in the Camden borough. These building geometries are part of the static input for the PALM model. Data for the shape and height of buildings are from Ordnance Survey (Ordnance Survey (GB), 2019b, a). Background map is OpenStreetMap (© OpenStreetMap Distributed under the Open Data Commons Open Database License (ODbL) v1.0 (for further information, please see OpenStreetMap)).

3.2 Anthropogenic emission datasets for London

For the London case study, anthropogenic CO₂ emissions are derived from the London Atmospheric Emissions Inventory (LAEI) (London Datastore, 2023), published by the Greater London Authority. LAEI employs a bottom-up methodology, combining activity data with sector-specific emission factors to estimate emissions across multiple sectors. The most recent dataset available during this study is LAEI2019, released in 2022, which reports annual total emissions for the year 2019 on a 1 km² grid covering the region within the M25 motorway. The inventory grid conforms to borough boundaries, so that grid cells intersecting multiple boroughs are subdivided accordingly.

Since LAEI only reports annual total emission on 1 km² tiles, it is essential to further disaggregate the emission data spatially and temporarily for use in an urban model. A direct approach would be to interpolate the LAEI 1 km² tile data on the finite difference grid defined in the urban model. However, this approach cannot capture the real emission scenario where emission within a neighbourhood can originated from several large localised emission sources such as gas boilers within households or industrial chimney stacks. While it is ideal to disaggregate all sectors into a fine grid map, there is a lack of available geo-spatial surrogate data for the task. For example, one might be able to find domestic building information in a city to disaggregate domestic energy usage sector, but there rarely exists a comprehensive list of construction site for non-road mobile machinery sector. Therefore, one needs to decide which sectors should be disaggregated with additional surrogate data and which should be simply interpolated. Here, we examine the proportions of each sub-subsector in the LAEI dataset. A list of the sector categories is included in Appendix A. Note that some sub-subsectors have been merged due to their small contribution. For sub-subsectors whose proportions are less than 10 %, their emissions would be simply distributed evenly in the finer grid and in time (annual total amount spread to hourly emission). Otherwise, we consider them to be significant and require fine-scale geospatial data as spatial proxies to disaggregate. In LAEI, CO₂ emission is reported for four major sectors (domestic, commercial and industrial, transport, and others) with sub-sectors and sub-subsectors. Figure 4 shows the significant emission sub-subsectors, which are domestic gas combustion, non-domestic (commercial and industrial) gas combustion, road transport, and industrial processes A1 (A1 refers to large industrial point sources).

Figure 4 Major CO₂ emission sub-subsectors in the Greater London Authority region in 2019 – annual emission of sub-subsectors and their proportion in total emissions.

It is noted here that although the word “disaggregation” is used to describe the operation of processing emission inventory data to the desired grid resolution (10 m), data in LAEI for the following sub-subsectors are not used: domestic gas combustion, non-domestic gas combustion, and industrial processes part A1. Instead, the source data for activities that LAEI has used to generate their 1 km² grid map are retrieved. Together with additional geospatial data as spatial proxies, we processed the 10 by 10 m² gird map for these sub-subsectors. Road transport was provided in LAEI as shapefile of the major road network and is convenient to interpolate to a grid of 10 by 10 m².

3.2.1 Spatial disaggregation

Spatial disaggregation of four significant sub-subsectors uses additional geospatial data. These include gas usage within Lower-layer Super Output Areas (LSOAs) or Middle-layer Super Output Areas (MSOAs), industrial point sources inventory, and domestic and non-domestic properties' energy performance certificates (EPC). Each LSOA is a statistical area covering 400 to 1200 households and has 1000 to 3000 residents, while each MSOA contains around 5000–15 000 residents or 2000–6000 households.

For domestic gas combustion, activity data comes from official statistics published by Department for Energy Security and Net Zero in the UK (Department for Energy Security and Net Zero, 2015). The data repository records gas usage (in kWh) in LSOAs. To disaggregate gas usage in LSOAs to households, we obtain property records from Energy Performance of Buildings Data (Department for Levelling Up, Housing & Communities, 2025). Each record in the dataset includes an estimate of annual CO₂ emissions for the property which are used as a spatial proxy. This estimation is calculated from many characteristics of a domestic property, for example, the floor area, the floor height, its main heating fuel, and the status of wall and window. We filter domestic properties that use gas as the main fuel, and then utilise their annual CO₂ emissions as a spatial proxy to distribute the total gas usage in an LSOA to properties. Then, the sum of gas usage from domestic properties that falls within a title of 10 m by 10 m is taken as the gas usage for the tile. The spatial join operation between the domestic property dataset and the gas usage dataset was performed using the PostgreSQL database with the PostGIS extension. The gas usage in tiles is multiplied by emission coefficient 0.20566 kgCO₂ kWh^{−1 2} to obtain annual emission in the grid map of 10 m by 10 m.

For non-domestic gas combustion, activity data also come from official statistics published by the Department for Energy Security and Net Zero in the UK (Department for Energy Security and Net Zero, 2015). Unlike domestic gas combustion, non-domestic gas usage statistics (measured in kWh) are reported at the MSOA level, which covers a larger geographical area than the LSOAs. To disaggregate gas usage in MSOAs to non-domestic properties, non-domestic energy records from Energy Performance of Buildings Data (Department for Levelling Up, Housing & Communities, 2025) are used, i.e. non-domestic EPC and DEC (display energy certificates). Unlike domestic records, many non-domestic records in this building dataset lack annual CO₂ emission estimation. Therefore, another floor_area attribute is used as a spatial proxy. Non-domestic properties that use gas as their main fuel are first filtered, and their floor areas are then used as a surrogate variable to apportion the total gas consumption within each MSOA to individual properties. The gas usage assigned to a 10 by 10 m² tile is calculated as the sum of the apportioned gas usage of all non-domestic properties intersecting that tile. The gas usage in tiles is multiplied by the emission coefficient 0.20455 kgCO₂ kWh^{−1 3} to obtain annual emission in the 10 by 10 m² grid-map. It is noted here that the non-domestic heating and power sub-subsector has excluded industrial energy usage in LAEI. Therefore, non-domestic gas combustion only contains those in commercial and industrial properties, i.e. heating and power in such buildings.

For road transport emissions, in addition to the 1 km² grid-map, LAEI presents a vector dataset of the London major road network stored in a shapefile. Each road section in the shapefile has an attribute for annual transport emission. The shapefile can be rasterised into a 10 by 10 m² emission grid-map using spatial joint operation. The residual emissions from road transport (which is called minor road transport emission in LAEI) represent about 10 % of the total road transport emission. These emissions are evenly distributed from the 1 km² grid to the 10 by 10 m² grid via interpolation.

The industrial processes part A1 in LAEI use data from both the Environment Agency's point source inventory (accessible from https://www.data.gov.uk/dataset/cfd94301-a2f2-48a2-9915-e477ca6d8b7e/pollution-inventory, last access: 15 July 2026) and the National Atmospheric Environmental Inventory (accessible from https://naei.energysecurity.gov.uk/data/maps/emissions-point-sources, last access: 15 July 2026). These inventories include power plants, wastewater processing factories, and other industrial sites. The point sources are spatially joined with the 10 by 10 m² tiles for a fine emission grid-map.

3.2.2 Temporal disaggregation

The temporal profiles of various sub-subsectors are used to disaggregate annual emissions into hourly emissions. Similarly as in spatial disaggregation, sectors with a proportion less than 1 % are disaggregated using constant temporal profiles. The other sectors use temporal profiles retrieved from CAMS-TEMPO (Guevara et al., 2021) and EDGAR (Crippa et al., 2020). These sectors and their temporal profiles are listed in Table 1.

Table 1 Sectors and their temporal profiles. In “Reference” column, CAMS refers to CAMS-TEMPO (Guevara et al., 2021).

Note that a different threshold (1 %) is selected for temporal disaggregation than in spatial disaggregation (10 %). The only additional sub-subsectors (between 1 % and 10 %) are: Industrial Oil/Coal Combustion, Industrial Construction NRMM, and Transport – Aviation. There are temporal profiles for these sub-subsectors in CAMS-TEMPO and EDGAR so we decide to include them in temporal disaggregation as well.

3.3 Vegetation indices and parameter selection

For the London case study, vegetation state variables required by the biogenic CO₂ flux model are derived from satellite-based surface reflectance products. Due to frequent cloud cover over the United Kingdom, EVI and LSWI are processed as monthly composites rather than at higher temporal resolution. For each calendar month in 2019, all available cloud-free satellite observations from the corresponding month throughout the 2018–2020 period are selected and temporally mosaicked. The median values are used to generate spatially complete monthly EVI and LSWI fields, reducing data gaps caused by cloud contamination.

The resulting monthly EVI and LSWI maps are provided to PALM as static input variables and are used by the biogenic flux module to represent vegetation phenology and water stress. This approach ensures spatial continuity of vegetation indices while preserving seasonal variability relevant for photosynthesis and respiration.

The parameters of the biogenic model are selected from the VPRM chemistry (WRF-Chem) weather research and forecast model (Skamarock et al., 2019), using parameter sets calibrated against measurements of the flux tower of eddy-covariance in Europe (Lian et al., 2021). Parameters are listed in Appendix F. These parameters are applied uniformly across the study domain, acknowledging that vegetation-specific calibration at the neighbourhood scale is not feasible given the available observations.

3.4 Diagnostic flux averaging configuration

To obtain simulated eddy-covariance (EC) fluxes from the LES, two approaches are employed. The primary approach computes EC fluxes from temporally averaged three-dimensional model fields using Eq. (10). An independent approach uses virtual measurement points that output high-frequency time-series data, from which EC fluxes are calculated following the standard eddy-covariance procedure. The agreement between these two approaches is used to verify the temporal-averaging method.

For the temporal-averaging approach, resolved vertical scalar fluxes are averaged over a 30 min window, consistent with the averaging period used in flux-tower observations. The mean vertical velocity and CO₂ concentration are obtained from the standard PALM output, while the covariance term is accumulated online during model integration. These correspond to the $\langle \stackrel{\mathrm{‾}}{w}\rangle$, $\langle \stackrel{\mathrm{‾}}{s}\rangle$, and $\langle \stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{s}\rangle$ terms in Eq. (10). All are 3D output quantities. These fields are then substituted into Eq. (10) in post-processing to obtain a 3D eddy-covariance flux field on all grid points. The 3D EC flux at the exact grid point where the BT tower measurement is located is extracted to compare the simulated EC flux with the site measurement described in the Sect. 3.5.

Besides the EC flux computed via Eq. (10), we also set up virtual measurement points in the domain to show the validity of this temporal averaging method. This is an independent calculation of the EC flux using time-series data. The virtual measurements output velocity component and CO₂ concentration at a frequency of 1 Hz, the same as that used in BT measurement (Helfter et al., 2016). After the simulation, 30 min means are computed from the 1 Hz time series and subtracted from instantaneous values to obtain the fluctuations w^′ and s^′. The EC flux is then calculated as the covariance $\langle w^{\prime }{s}^{\prime }\rangle$ over each half-hour period. The comparison of the EC fluxes obtained via virtual measurement and via temporal averaging (Eq. 10) is shown in Fig. D1. The two approaches yield almost the same flux profiles. In the main text, we use the EC flux computed from temporal averaging 3D fields as the simulated flux to compare with BT tower measurement.

A 30 min averaging window is used throughout this study to match the averaging period of the BT measurements and standard eddy-covariance processing. While EC fluxes are computed using 30 min averaging intervals, a 2 h moving average is applied to both simulated and measured flux time series when calculating RMSE and correlation coefficients to reduce high-frequency variability. The moving average is computed by taking four half-hourly EC fluxes and compute their mean. The timestamp for the moving average is at the end of the 2 h window.

3.5 BT Tower Eddy-Covariance Measurements and Source Area Characteristics

The measurement site is in Camden, London. An eddy-covariance flux tower was set up on the BT tower in Camden by previous researchers (Wood et al., 2010; Helfter et al., 2011). It collected the vertical CO₂ flux (Helfter et al., 2011), turbulent flow quantities (Wood et al., 2010) and other meteorological and greenhouse gas variables. We used the measurements in 2019 to validate our CO₂ model. Specifically, the wind speed and the eddy-covariance CO₂ flux were used in the validation. The temporal interval of measurements was half an hour. The CO₂ flux is computed via the procedure described in Sect. 2.5. In the following, we provide some essential information on the measurement site and its surrounding areas. For more details, the reader is referred to Helfter et al. (2016).

The eddy-covariance (EC) measurements of CO₂ were conducted from the BT Tower in central London ($\mathrm{51}\mathit{°}{\mathrm{31}}^{\prime }{\mathrm{17.4}}^{\prime \prime }$ N, $\mathrm{0}\mathit{°}{\mathrm{08}}^{\prime }{\mathrm{20.0}}^{\prime \prime }$ W). The EC system was installed on a lattice structure mounted on the roof of the telecommunications tower, providing an effective measurement height of approximately 190 m above street level. The system consisted of a three-dimensional sonic anemometer and a closed-path cavity ring-down spectrometer measuring CO₂, CH₄ and H₂O, with air sampled through a 45 m inlet line. Owing to the unusually large measurement height, the footprint of the observations extended over several kilometres and integrated emissions from a substantial portion of central London.

The BT Tower is located within a densely urbanized area of central London characterized by a mixture of commercial, residential and transport infrastructure. Within a radius of approximately 10 km, the mean building height is 8.8 ± 3.0 m, while typical suburban building heights are 5.6 ± 1.8 m. Consequently, the 190 m measurement height is more than an order of magnitude greater than the average surrounding building height, placing the EC sensors well above the urban roughness sublayer and enabling observation of integrated fluxes from a large urban source area.

Footprint analyses indicate that the source area varies seasonally from a few kilometres under convective conditions to several tens of kilometres during wintertime stable conditions. The footprint encompasses a diverse range of land uses. To the southwest and northwest, the source area includes two major urban parks, Hyde Park (142 ha) and Regent’s Park (197 ha), respectively. Northern sectors are dominated by suburban residential neighbourhoods, whereas eastern and southern sectors contain a mixture of densely built residential and commercial districts. The southeastern footprint also includes part of the River Thames. As a result, the measured CO₂ fluxes represent the integrated influence of emissions from road traffic, residential and commercial energy use, and human activities across a heterogeneous urban landscape.

3.6 Evaluation of biogenic module

Because CO₂ fluxes measured at the BT tower are dominated by anthropogenic emissions, they provide limited information on the realism of the biogenic component of the model. Therefore, an additional qualitative evaluations of the implemented VPRM scheme were made.

First, we conducted using eddy-covariance measurements from the ICOS site CZ-Lnz, a deciduous forest ecosystem in the Czech Republic (ICOS RI et al., 2023). The objective of this comparison is not to perform a site-specific validation of vegetation fluxes in Regent's Park, but rather to assess whether the model reproduces realistic magnitudes, seasonal variations, and diurnal patterns of net ecosystem exchange (NEE) for deciduous vegetation. Given the climatic and ecological differences between the two sites, the comparison is interpreted as an assessment of the VPRM implementation rather than as a direct validation of local biogenic fluxes.

The CZ-Lnz site is located in a deciduous floodplain forest, and its eddy-covariance measurements are primarily influenced by local biogenic carbon exchange. Although the climatic and ecological conditions differ from those in London, CZ-Lnz provides a long-term, high-quality record of net ecosystem exchange (NEE) from a deciduous forest ecosystem, allowing a qualitative assessment of the VPRM representation of forest carbon uptake and respiration.

Monthly ensemble averages of tower observations were used to derive representative half-hourly diurnal NEE profiles for each month. For the model, NEE was extracted directly from the VPRM output for a grid cell in Regent's Park classified as deciduous forest. The comparison therefore focuses on the temporal characteristics of NEE rather than on direct site-to-site agreement.

Second, the total NEE for the parent domain (roughly Greater London Authority area) was compared with the ICOS-VPRM product (Gerbig and Koch, 2021). This product runs offline VPRM module with MODIS-derived satelite indices and meteorological inputs from ECMWF IFS analysis data. The product has a spatial resolution of ∼ 0.1°. The daily sum of NEE for the twelve simulated days from PALM-CO₂ and from ICOS-VPRM were compared. Using the simulated days as representative values for their corresponding months, annual NEE rate (both total amount and monthly correlation) were reported.

3.7 Modelling settings in PALM-CO₂

3.7.1 Static driver

As an urban model, PALM takes all the data describing urban land surface within one file called the static driver. These data include terrain height, land-use type and building geometries. The data sources for static inputs in the study area are listed in Table 2. The static input file is prepared with an open source tool GEO4PALM (Lin et al., 2024).

For buildings' geometry, the topography map from Ordnance Survey (OS) (Ordnance Survey (GB), 2019b) is used. The full topography map is filtered to include buildings only. For buildings' height, an additional OS building height dataset (Ordnance Survey (GB), 2019a) is used. The shapes (geometries) and the heights of buildings are matched with unique IDs indexed by OS to form a vectorised building map. For terrain height, two data sources of different resolution are selected for the parent and child domain, namely OS Terrain 50 (Ordnance Survey (GB), 2019c) (resolution at 50 m) and National Lidar Programme (Environment Agency, 2025) (resolution at 1 m). Land-use and vegetation types are retrieved from ESA WorldCover dataset (Zanaga et al., 2022).

Using these static input data, we activated the following surface modules. An urban surface module solves for the buildings' energy balance and energy exchange with outdoor environments. A land surface module computes the drag and heat dispersion in different land-use types, including vegetation, pavement, water, and soil. Land-use types also affect radiation from and into the urban surface. The soil temperature is modelled with an eight-layer energy balance equation.

3.7.2 Initial and boundary conditions (dynamic driver)

The initial and boundary conditions are taken from real meteorological data from the ECMWF reanalysis dataset (ERA5) (Hersbach et al., 2020). These are prepared as a “dynamic driver” file in PALM. ERA5 model level data for hourly velocity, potential temperature, humidity, long- and short-wave radiation, and surface level pressure are used to drive the urban model simulation.

For CO₂ transport, the background value of its mixing ratio is also required on the boundary outside the modelling region. This information comes from CAMS global greenhouse gas reanalysis (EGG4) (Copernicus Atmosphere Monitoring Service, 2021), due to its spatially and temporally consistent global coverage and suitability for boundary forcing in regional atmospheric transport simulations. 3-hourly CO₂ mixing ratio in the study region is used as the background value on the boundary. It is interpolated into hourly data and written in the dynamic file. These dynamic driven data sources are also listed in Table 2.

(Ordnance Survey (GB), 2019c)(Environment Agency, 2025)(Ordnance Survey (GB), 2019b)(Ordnance Survey (GB), 2019a)(Zanaga et al., 2022)(European Space Agency(ESA), 2025)(Hersbach et al., 2020)(Copernicus Atmosphere Monitoring Service, 2021)

Table 2 Data sources for static and dynamic drivers in the study area.

3.7.3 Running on HPC

The case is set up to running on ARCHER2 (UK HPC) using 1274 cores, out of which 250 cores are used for the parent domain (N01) with 10 by 25 subdomains in the x and y directions and the rest 1024 cores are used for the child domain (N02) with 32 by 32 subdomains in the x and y directions. The 16th day of each month in 2019 is selected for the simulation at the study area. Depending on the average wind speed of the day, the CPU wall time for one day simulation takes between 6.1 h (for 16 April) and 23.2 h (for 16 March).

4 Results and discussion

4.1 Validation against BT tower measurement

In this subsection, we show comparison of the modelled flux and the measured flux by the eddy covariance method. This comparison validates: (1) the input anthropogenic emissions and online-computed biogenic fluxes, and (2) the CO₂ transport process calculated by PALM.

In previous studies, CO₂ transport has been approximated using simplified two-dimensional advection–dispersion approaches, in which turbulent source areas or flux footprints were parameterized rather than explicitly resolved; see Christen et al. (2011) for an urban carbon emission model that leverages the 2D model for validation. In contrast, PALM resolves 3D turbulence using large-eddy simulation, enabling a more detailed representation of turbulent dispersion and vertical transportation processes, and thereby providing a more physically consistent characterization of source-receptor relationships.

To reduce the influence of high-frequency variability and facilitate comparison between simulated and measured fluxes, both time series were additionally smoothed using a 2 h moving average prior to the calculation of RMSE and correlation coefficients. The 30 min EC fluxes were retained for all other qualitative analyses. Wind speed and wind direction have also been smoothed using the same method.

Figure 5 shows the validation results for the CO₂ flux for selected months representing each season of a year. Overall, good agreement between the simulated CO₂ flux and the measurement at the BT tower has been reached. When July is excluded, the average RMSE and correlation coefficient are 16.74 µmol m⁻² s⁻¹ and 0.74, respectively. The trend of the diurnal cycle has been captured by the model and matches well with observations. Specifically, the CO₂ flux remains low in the early morning before sunrise, increases during the morning hours (08:00 a.m. to 12:00 p.m.), reaches a peak in the afternoon (04:00 p.m. to 06:00 p.m.), and gradually decreases toward midnight (00:00 a.m.). As the measurement site is located in a commercial-residential mixed area with limited green spaces, the flux pattern is dominated by anthropogenic emission. Therefore, such a diurnal cycle meets the expectation of human activity patterns within a day. Quantitatively, the magnitudes of fluxes in January, April, and October show good agreement with observations. Statistics show relatively low error and high correlation between simulation and measurement.

Figure 5 Comparison of modelled CO₂ fluxes and BT tower measurements in 2019. Dotted lines and markers show original temporal data at half-hourly resolution for both PALM model and tower measurements. Solid lines show two-hour moving average. All times in the figure are in local time. Sporadic missing data in measurement have been filled with interpolated values. (a) 16 January 2019; (b) 16 April 2019; (c) 16 July 2019; (d) 16 October 2019.

Modelled fluxes in July show a consistent pattern as in other seasons, while discrepancy from measurement exists. The site measurement in 16 July has poor quality which already shows a non-typical diurnal pattern in the figure. Besides, an overestimation of flux on 16 July persists. As one would expect that anthropogenic emission in July (summer time) is smaller than in October and January, yet the simulated flux in July is higher. This overestimation originates from extremely weakly forced atmosphere during daytime peak hours. As can be seen in Fig. 6c, wind speed at 190 m height is close to 1 m s⁻¹ between 10:00 to 16:00. Poor advection caused build-up of CO₂ in the domain and an over-estimation of EC flux. The quality control flag and additional analysis of July data can be found in Appendix C.

Figure 6 Comparison of modelled wind speed and BT tower measurements in 2019. Dotted lines and markers show half-hourly averaged wind speed from simulation and measurement. Solid lines show two-hour moving average. Sporadic missing data in measurement have been filled with interpolated values. (a) 16 January 2019; (b) 16 April 2019; (c) 16 July 2019; (d) 16 October 2019.

Figures 6 and 7 show validation of the wind speed and direction. The average RMSEs for wind speed and wind direction in the figures were 2.30 m s⁻¹ and 37.50°, respectively, while the corresponding average correlation coefficients were 0.50 and 0.68. The wind velocity on the 10 by 10 m² grid at 190 m height corresponding to the measurement location is extracted from the simulation. In the figures, wind speed and direction at a half-hourly interval and its 2 h moving average daily profile are shown. The general trend of daily variation of wind (the 2 h moving average profile) matches that in the measurement. For all days simulated in 2019, statistics show good agreement between simulation and measurement for both wind speed and direction, with average RSME around 2.4 m s⁻¹ and 28.65° and correlation coefficient 0.55 and 0.75 (see Appendix B). The validation results demonstrate that ERA5 provides suitable dynamic boundary conditions for the parent domain and confirm the effectiveness of the nesting strategy in representing processes within the nested domain.

Figure 7 Comparison of modelled wind direction and BT tower measurements in 2019. Dotted lines and markers show half-hourly averaged wind speed from simulation and measurement. Solid lines show two-hour moving average. Sporadic missing data in measurement have been filled with interpolated values. (a) 16 January 2019; (b) 16 April 2019; (c) 16 July 2019; (d) 16 October 2019.

Results for all months can be found in Appendix B.

4.2 Evaluation of biogenic fluxes

Comparison with the BT tower measurements provides the primary validation of the modelling framework in the study area. However, CO₂ fluxes measured at the BT tower are dominated by anthropogenic emissions, making it difficult to directly assess the realism of the vegetation component of the model. Here, we compare biogenic flux with two independent sources in the literature to evaluate the biogenic fluxes in PALM-CO₂.

To provide an independent evaluation of the implemented VPRM scheme, we first compare the simulated biogenic fluxes against eddy-covariance observations from a European forest site, Lanzhot in Czech (CZ-Lnz) (ICOS RI et al., 2023). The purpose of this comparison is not to perform a site-specific validation of vegetation in Regent's Park, but rather to assess whether the model produces realistic magnitudes, seasonal variations, and diurnal cycles of ecosystem carbon exchange.

Figure 8 compares the simulated and observed NEE profile. Both datasets exhibit the expected seasonal behaviour of temperate deciduous forests. During winter (December, January, February), NEE remains relatively small and positive, reflecting limited photosynthetic activity and continued ecosystem respiration. During summer (May, June, July), stronger daytime carbon uptake produces larger negative NEE values and a pronounced diurnal cycle. The model reproduces these broad seasonal and diurnal characteristics, indicating that the implemented VPRM module captures the dominant controls on biogenic carbon exchange.

Figure 8 Comparison of monthly 24 h net ecosystem exchange profile from flux tower measurement and from PALM model. The vegetation type is deciduous trees. The measurement location is at Lanžhot in Czech Republic (CZ-Lnz). Shaded area is the standard deviation range for measurement. Vertical dashed lines denote month boundaries.

Differences are evident during the spring green-up period and autumn senescence. The simulated NEE shows stronger carbon uptake in March and April and weaker uptake in August, September and October compared with CZ-Lnz observations. These discrepancies are likely related to differences in phenology between the two sites since they correspond to the budburst and senescence stages of deciduous forests (Duchemin et al., 1999). The timing of budburst and senescence is influenced by local climate, species composition, soil conditions, and management practices, all of which differ between the urban park environment of Regent's Park and the rural floodplain forest represented by CZ-Lnz. As shown in Fig. 9, EVI and LSWI increase rapidly in March at the Regent's Park location, leading to enhanced simulated GPP through Eq. (2).

Figure 9 Comparison of monthly EVI and LSWI at two deciduous forests site: Regent Park in Camden, and Lanžhot in Czech Republic (CZ-Lnz). Biogenic fluxes in both sites are compred in Figure 8.

The stronger uptake simulated during spring is consistent with the rapid increase in EVI and LSWI at the Regent's Park location. However, the differences observed during late summer and early autumn cannot be explained by EVI alone. Although EVI remains relatively high during August and September, NEE is influenced by additional environmental controls in VPRM, including temperature and moisture stress, as well as site-specific ecosystem characteristics. Furthermore, the comparison involves two ecologically distinct systems: an urban park in London and a rural floodplain forest in the Czech Republic. Differences in species composition, soil moisture conditions, forest structure, and phenological timing may therefore lead to different seasonal trajectories of carbon uptake despite similar vegetation index values. Consequently, the late-summer discrepancy is likely attributable to site-specific ecological differences rather than changes in canopy greenness alone.

Second, the sum of the total biogenic flux in the N01 domain (60 by 60 km²) is compared with the ICOS-VPRM product (Gerbig and Koch, 2021). Figure 10 shows the comparison. The annual trend of biogenic flux computed in PALM agrees with that in ICOS-VPRM in general. A positive net flux during winter, late summer and autumn seasons appeared in both computations, while spring and early summer show negative net flux. The values are listed in Table 3. Using the NEE on the 16th day of each month to represent monthly flux, average biogenic flux in 2019 is −3.95 kton d⁻¹ in our model while ICOS-VPRM reports average NEE −4.92 kton d⁻¹. Correlation coefficient between them is 0.97. To show the ratio of biogenic fluxes and anthropogenic emissions, Table 3 also lists the total biogenic fluxes in the N02 domain (8 by 8 km²) as well as the N01 domain. For the N02 domain, biogenic flux to anthropogenic flux ratio is between −4.72 % to 2.29 %; for the N01 domain the ratio is between −78.2 % and 8.6 %. This shows that in the densely populated Camden borough biogenic fluxes have very limited offset to anthropogenic emissions.

Figure 10 Comparison of daily total biogenic flux (NEE) computed by PALM and in the ICOS-VPRM product. Total fluxes in the N01 domain (60 by 60 km²) is plotted. For ICOS-VPRM, the area is approximated as the grid does not align exactly with the BNG projection used in PALM.

Table 3 Daily integrated NEE and anthropogenic CO₂ emissions on the 16th day of each month in 2019. NEE values are reported for the PALM N02 domain (8 × 8 km²), the PALM N01 domain (60 × 60 km²), and the ICOS-VPRM product over approximately the N01 domain.

4.3 Urban boundary flow and carbon dioxide concentration

In this subsection, we show results on carbon dioxide concentration and its relation to the height of the Atmospheric Boundary Layer (ABL). Three snapshots of CO₂ concentration in the morning 08:30 a.m., noon 12:30 a.m., and evening 08:30 p.m. are selected in the following analysis. The ABL heights and Obukhov lengths for the snapshots are listed in Table 4. Figure 11 shows the CO₂ concentration at the surface level at selected times during a day in January and July. Due to terrain height and building geometries, the surface level at different horizontal places are at different levels. This has been processed and in the figure only the lowest level in the air is plotted. In all three snapshots in both January and July, individual large anthropogenic emission locations such as major roads and buildings in the commercial region (lower right in each sub-figures) can be identified. They contribute substantially to the variation of local concentrations. In January, due to a relatively higher wind speed (about 10 m s⁻¹, see Fig. 6a), advection near the ground level is prominent throughout the day. However, its effect weakens around noon (12:30 p.m.), due to an increase in upward turbulent flow driven by ground heating in the air. At night, reduced heat turbulence causes a build-up of CO₂ in the region. In July, the wind speed is lower (about 4 m s⁻¹, see Fig. 6c). In the morning, only limited advection occurs due to weak lateral wind and minimal heat turbulence. Therefore, a large concentration region appears in the dense commercial region (lower right in the sub-figure). At noon, due to increased heat turbulence and upward flow, concentration is reduced in the region. At night, mixed effects from suppressed upward turbulence and increased lateral advection (due to increased wind speed) produce a moderate concentration. 24 h surface concentration profiles for the two days discussed here can be found in Appendix G.

Figure 11 Instantaneous surface level CO₂ concentration at different times in January (a–c) and July (d–f) from the LES simulation.

Table 4 ABL height statistics and Obukhov length for 16 January 2019 and 16 July 2019.

The January cases are characterized by a relatively shallow boundary layer (411–827 m) and near-neutral to weakly unstable/stable conditions ($\left|L\right|\approx$ 100–300 m). In contrast, the July cases exhibit substantially stronger instability during daytime (L=−1.9 to −7.4 m), leading to enhanced turbulent mixing and the development of a deep convective boundary layer reaching 1.7 km at midday. Following sunset, both periods transition to stable stratification (L>0), accompanied by a reduction in boundary-layer height, although the summer boundary layer remains considerably deeper than the winter boundary layer due to stronger daytime heating and residual turbulence. The ABL height and stability are reflected in the vertical CO₂ concentration profile. Figure 12 shows the concentration on a x=4000 m vertical cross-section over the domain at different times in January and July. The location of the vertical cross-section is shown in Fig. 2b. In general, during morning and night, the shallow boundary layer traps CO₂ below approximately ABL height, while during the day the ABL is thicker and the unstable ABL can mix CO₂ well. An exception occurs in the morning (08:30 a.m.) in January when advective transport lifts some plumes above the boundary layer, reducing near-surface accumulation.

Figure 12 Instantaneous CO₂ concentration at x=4000 m cross-section in the domain at different times in January (a–c) and July (d–f). LES results are shown here. On the bottom the annotated axis illustrates the topography and geographic references along the cross-section. For ABL heights and stability corresponding to the snapshots see Table 4. The location of the vertical cross-section is shown in Fig. 2b.

Figure 12 also shows the impact of spatial patterns of anthropogenic emissions on CO₂ concentrations within the lower atmospheric layer. As shown in the geography tags, urban land-use types are interwoven from south to north (+Y direction in the figure), alternating commercial districts and green spaces: beginning with parkland (Hyde Park), followed by dense commercial areas (Marylebone), then another parkland and commercial areas (Regent's Park and Camden Town), mixed natural land and residential zones (Hampstead Heath and Highgate). Concentration hotspots correspond to densely populated areas with higher human activities. Throughout the day and in both months, the effect of turbulent dispersion (plumes) can be clearly observed in the figure.

The above pattern is revealed through LES that explicitly represents urban topography and spatially heterogeneous anthropogenic and biogenic fluxes at a 10 m resolution. The resulting complex CO₂ dispersion pattern highlights the importance of resolving urban-scale flow structures and emission heterogeneity when simulating urban carbon transport. On this basis, the model can help constrain emission patterns in an inversion problem. This conclusion echoes the argument in Brunner et al. (2019) that vertical distribution of carbon fluxes greatly affects the transport process.

4.4 Grid sensitivity study

The grid resolution in the case study was selected as 10 m in the child domain. This is a balance between resolution, availability of static and dynamic driver data, and the purpose of the carbon dioxide case study in general. From the perspective of numerical discretisation, one wishes to have a high resolution to resolve turbulent transport at as small scales as possible. However, this can be prohibitive for computational cost. It also suggests a new requirement for static and dynamic drivers. For static driver, it requires finer topography data (buildings and terrain height), land use data, and vegetation data which could be scarce. Yet from the perspective of the user end, a study on carbon dioxide does not require as fine a resolution as those used in air pollution study where the individual exposure can benefit from a finer concentration map. Eventually, we selected the 10 m resolution in the main case study.

Nevertheless, it is worth exploring the turbulent transport process at different grid resolution. Here, we include two additional simulations of an ideal convective boundary layer (CBL) case at 10 and 5 m resolution and compare the results. The domain is 2 km by 2 km by 3.5 km, i.e. a smaller area surrounding the BT tower in the child domain in the main case study shown in Fig. 2. The 2 km domain is created using the same Ordnance Survey data for terrain height and buildings at 10 and 5 m resolution while the flow and emission conditions are controlled. The flow is driven by both the geostrophic wind, u_g=5 m s⁻¹ and v_g=5 m s⁻¹ and the bottom heating at 0.1 K m s⁻¹. These values are typical in the study area. The side boundaries for momentum and temperature are cyclic. The top boundary for momentum is the free slip, while for energy it is a constant gradient at 0.01 K m⁻¹. The emission of passive scalar CO₂ is set as a constant rate in time but spatially heterogeneous, whose pattern is determined by a typical winter day in the study area. The side boundaries for CO₂ are set as 0, while the top boundary is zero gradient Neumann. The initial flow field is geostrophic wind. The initial potential temperature is constant below 800 m and increasing at 0.01 K m⁻¹ until the top boundary to form an inversion layer. The initial concentration for CO₂ is set as 0. Both cases have been run for 48 h. Horizontally, the grid resolution is 10 and 5 m respectively. Vertically, below 400 m the grid resolution is 10 and 5 m respectively, while above 400 m the grid space increases at a ratio of 1.08 until reaches approximately 3.5 km.

Figure 13 shows the vertical flux profile at the end of simulation (48 h) for two grid sizes. Note that the temporal mean is used to compute the turbulent fluctuation and turbulent flux for subplot (a) and (e), while spatial (horizontal) mean is used for the purpose in the rest subplots. In Fig. 13b and f, resolved fluctuations are defined using horizontal averaging. As a result, the resolved flux may include contributions from persistent spatial heterogeneity arising from terrain-induced flow structures and spatially varying emission patterns. Therefore, the resolved flux fraction should be interpreted as the portion of transport resolved by the LES grid, rather than as a strict decomposition between turbulent and dispersive transport processes. This interpretation becomes more reliable above the roughness sublayer, where turbulence is increasingly horizontally homogeneous and the influence of dispersive motions associated with surface heterogeneity is diminished. In the current case, average building height or roughness height is about 37.7 m for 10 m grid resolution and 33.2 m for 5 m grid resolution.

Overall, the flux profiles obtained with the two grid resolutions show good agreement at intermediate and higher elevations, particularly above approximately 75 m (about two roughness heights). Below this level, however, notable differences emerge. As shown in Fig. 13a and e, the turbulent flux at 5 m resolution exhibits a larger standard deviation, reflecting the stronger influence of surface heterogeneity captured by the finer representation of urban morphology, including street canyons, building geometry, and spatially varying emission patterns.

Figure 13 Comparison of vertical turbulent flux profile at 48 h of the 10 and 5 m resolution simulations. (a) The horizontal mean (solid line) and standard deviation range (shaded area) of turbulent flux computed from temporal average Eq. (12) on each vertical level; (b) resolved turbulent flux computed from spatial average; (c) under-resolved turbulent flux computed from spatial average; (d) percentage of resolved flux. (e)–(f): zoomed-in view of (a)–(d) for 0 to 200 m height.

The extent to which turbulent transport is resolved also differs within the roughness sublayer, as illustrated by the profiles of resolved and under-resolved turbulent fluxes in panels b–d and f–g. Below 75 m, the 5 m simulation resolves a greater fraction of the turbulent flux than the 10 m simulation. At higher elevations (> 100 m, or roughly three roughness heights), both resolutions resolve the vast majority of the turbulent transport, with resolved flux fractions exceeding 96.9 % and 98.5 % for the 10 and 5 m simulations, respectively. At 190 m, corresponding to the height of the BT Tower measurements, the resolved flux fractions are 98.9 % for both grid resolutions.

These results indicate that, under convective boundary-layer conditions, a 5 m grid resolution provides a more accurate representation of surface heterogeneity and turbulent transport within the lower urban canopy and roughness sublayer (up to approximately three roughness heights). Above this region, however, turbulent transport of passive scalars is already well resolved at a 10 m grid resolution. Additional details of the grid-sensitivity experiment and model configuration are provided in the Appendix D.

4.5 Ensemble simulation of a representative validation case

While the validation in previous subsections show promising agreement, it should be noted that the LES presented there represents only one realisation of turbulent flow field for a given day. Because atmospheric turbulence is chaotic, even under identical large-scale forcing and boundary conditions, small perturbations in the turbulent field can lead to different instantaneous evolutions. Therefore, a perfect match between a single LES realization and the measured half-hourly eddy-covariance fluxes should not be expected. To better quantify the variability associated with turbulent realization uncertainty, a 10-member ensemble simulation was performed for 16 October 2019.

The ensemble members were generated by using different random number seeds in the PALM synthetic turbulence generator, while keeping all other model settings unchanged. In this way, the ensemble isolates the sensitivity of the simulated EC flux to differences in the turbulent realization, rather than to changes in forcing, domain configuration, or model physics.

Figure 14 shows the half-hourly EC flux from all 10 ensemble members together with the ensemble mean. At the native 30 min resolution, both the model output and the measurements exhibit strong short-term fluctuations. Such variability is expected for turbulent fluxes and makes a direct point-by-point comparison between simulation and observation difficult, especially when the purpose is to assess whether the LES reproduces the overall temporal evolution and magnitude of the fluxes. For this reason, the comparison was additionally carried out using a 2 h moving average, applied consistently to both the LES output and the measurements. This smoothing reduces the influence of high-frequency fluctuations while retaining the diurnal structure of the flux evolution.

Figure 14 Ensemble members for 16 October 2019 and their mean. Top plot shows original half-hourly EC flux processed from LES; bottom plot shows 2 h moving averaged EC flux. Thin coloured lines represent individual ensemble members. Thick black lines represent their means.

The smoothed ensemble results indicate that, despite substantial spread among individual members at half-hourly resolution, the realizations show a broadly consistent diurnal pattern after applying the 2 h moving average. The ensemble mean captures the main daytime increase in EC flux, with values rising from low nighttime and early-morning levels to a midday/afternoon maximum, followed by a decline toward evening. This suggests that the large-scale temporal development of the flux is robust across different turbulent realisations, whereas part of the short-term mismatch between a single simulation and the observations can be attributed to intrinsic flow variability.

Figure 15 compares the measured EC flux and the ensemble mean using the 2 h moving average. The black line represents the ensemble mean, and the shaded band indicates the uncertainty range of the ensemble mean given by the standard deviation. The red line shows the smoothed measurements. Overall, the ensemble mean reproduces the broad diurnal evolution seen in the observations, including the increase during the morning, elevated daytime fluxes, and the decrease in the evening. Quantitatively, the comparison yields an RMSE of 16.60 µmol m⁻² s⁻¹, and a Pearson correlation coefficient of 0.77 for the 2 h moving-mean series, indicating moderate-to-good agreement in the temporal evolution.

Figure 15 Comparison of ensemble mean and BT tower measurement on 16 October 2019. Shaded area shows the standard deviation range of the ensemble.

At the same time, differences remain between the modelled and observed fluxes, particularly during periods when the measured half-hourly values show sharp excursions. Some of these discrepancies are reduced after temporal smoothing, which supports the interpretation that part of the disagreement in the single-realisation comparison arises from the stochastic nature of turbulence rather than from a systematic model error alone. However, the remaining differences also indicate that ensemble variability does not explain all mismatches between LES and observations. The residual mismatches can be due to uncertainties in the emission profiles of the particular day.

5 Conclusions

We developed PALM-CO₂, a high-resolution urban carbon dioxide transport modelling framework based on an urban large-eddy simulation code, PALM. The model resolves CO₂ transport at metre-scale resolution (down to 10 m) under realistic meteorological conditions and explicitly represents urban topography, including terrain, buildings, and vegetation types. A key novelty of PALM-CO₂, is its fully coupled treatment of fine-scale atmospheric dynamics with spatially detailed anthropogenic and biogenic surface carbon fluxes, enabling the simulation of CO₂ transport processes within the urban canopy layer that cannot be captured by coarser-scale models.

The model considers complex urban land surface carbon fluxes from both anthropogenic and biogenic activities. Anthropogenic CO₂ emissions are prescribed at metre-scale resolution (i.e. 10 m) using a bottom-up disaggregation framework that integrates multiple urban data sources, including gas consumption statistics, building energy performance data, industrial point-source inventories, and road transport emissions. Biogenic CO₂ fluxes are represented through the online implementation of the VPRM within PALM. VPRM dynamically computes net ecosystem exchange using PALM-resolved near-surface meteorological conditions and satellite-derived vegetation indices (EVI and LSWI), allowing for a physically consistent coupling between surface carbon exchange and resolved turbulence.

The model is validated against eddy-covariance flux and wind measurements from measurement site at BT tower. Within the urban case study domain, simulated turbulent CO₂ fluxes and horizontal wind speeds show good agreement with observations from the BT Tower, where anthropogenic emissions dominate the measured signal. When July is excluded, the average RMSE and correlation coefficient for EC fluxes are 16.74 µmol m⁻² s⁻¹ and 0.74, respectively. The average RMSEs for wind speed and wind direction were 2.30 m s⁻¹ and 37.50°, respectively, while the corresponding average correlation coefficients were 0.50 and 0.68. To further assess the biogenic component, PALM-CO₂ is evaluated against flux measurements from the Lanžhot site in the Czech Republic, a deciduous forest ecosystem comparable to the green spaces in the Regent's Park in the study area. Simulated NEE captures both dormant and peak growing season behaviour well, while differences in the timing of budburst and senescence reflect climatic contrasts between the two regions, as also indicated by EVI and LSWI. The annual average biogenic flux of the Greater London Area is −3.95 kton d⁻¹ in the model, close to the number in the ICOS-VPRM product −4.92 kton d⁻¹. Annual cycle of biogenic flux agrees with the ICOS-VPRM product, with a correlation coefficient 0.97. These evaluations demonstrate the robustness of both the transport and surface flux components of PALM-CO₂.

Overall, PALM-CO₂ provides a physically consistent framework for simulating carbon dioxide transport in complex urban environments at metre-scale resolution. By explicitly resolving urban topography, atmospheric turbulence, and the coupled effects of anthropogenic and biogenic surface fluxes, the model captures multi-scale interactions that cannot be represented in coarser atmospheric transport models. The results demonstrate that fine-scale heterogeneity in emissions, land cover, and boundary-layer dynamics strongly influences urban CO₂ distributions, particularly under stable conditions and during periods of shallow boundary-layer development. These findings underline the importance of high-resolution large-eddy simulation approaches for advancing process-level understanding of urban carbon dynamics and for improving the reliability of urban-scale carbon modelling. PALM-CO₂ thus establishes a robust modelling foundation for future studies investigating urban carbon processes, emission attribution, and the role of cities in the regional carbon cycle.

Appendix A: Sectors and subsectors in emission inventory

Table A1 lists the sectors, sub-sectors and sub-subsectors used in the London Atmospheric Emission Inventory (London Datastore, 2023). Sub-sectors or sub-subsectors with no emissions have been excluded in the table. Some sub-subsectors have not shown in the table and merged (those of Rail and River).

Table A1 Inventory source categories and GLA total emissions. NRMM = non-road mobile machinery.

Appendix B: Results for all months

Results for all months are presented in Figs. B1, B2, and B3 for EC fluxes, wind speed, and wind direction respectively. The metrics are listed in Table B1.

Table B1 Monthly performance metrics for EC flux, wind speed, and wind direction predictions. Root mean square error and correlation coefficient r are listed. Units are for RSMEs and biases. r has unit 1. Both simulation and BT tower measurement are on the 16th day of each month. Measurement is missing for 16 June.

Figure B1 Comparison of simulated and measured EC fluxes for all days. Sporadic missing values in BT tower measurement have been filled by linear interpolation. All day's data are missing for 16 June 2019.

Figure B2 Comparison of simulated and measured wind magnitude.

Figure B3 Comparison of simulated and measured wind direction.

Appendix C: Additional analysis for 16 July

LES shows large over-estimate of EC flux on 16 July in Fig. 5. The reason is two-fold. First, the measurement on 16 July is not representative of a typical day in July. The processed EC flux data in the BT tower dataset have gone through quality control. The quality control checks several criteria including the magnitude of turbulence intensity, number of valid instantaneous measurements, out-of-range velocity ratios and range of heat fluxes. A quality control flag (qc) is provided in the dataset (Helfter et al., 2016). A qc value less than 2 is deemed good data. Figure C1 shows the BT measurement EC flux in July 2019 and their individual qc flag. It can be seen that from 16 July, quality of data deteriorated and more data are with qc flag = 2. The plot also shows the average 24 h EC flux profile for the first two weeks of July (1 to 14) as a blue line in the first day. This average profile is very different to that on 16 July. Therefore, it is likely the measurements on 16 July are unfortunately not representative of the real situation.

Another reason for the over-estimation in July is the small wind speed during morning peak hours on 16 July. As can be seen in Fig. 6 in main text, wind speed from 08:00 to 18:00 is very low even at 190 m height. At lower ground the wind speed would be even smaller. The worse convection condition will likely cause the build-up of flux in the domain that would otherwise be transported horizontally out of the domain. Eventually this caused an over-estimation.

To further prove this point, we select another day in July (the 20th) with moderate wind magnitude. This day has similar anthropogenic emission to 16 July. The results are plotted in Fig. C2. The simulated EC flux (blue lines in figure) is smaller, and smaller than that in January and October as expected. Comparison with observation is unreliable for 20 July due to bad data quality (orange lines in figure). However, if the average 24 h profile for the first two weeks of July (green lines in figure) is used in the comparison, the 20 July simulation is in good agreement with measurements (RMSE = 8.28 umol m⁻² s⁻¹, r= 0.819).

Figure C1 BT measurement and quality control flag for July 2019. Individual half-hourly EC flux data points are coloured with quality flag. Less than 2 means good quality.

Figure C2 Comparison of EC flux from LES and from measurement on 20 July. Due to poor quality of measurement on 20 July, the average EC flux profile of 1 to 14 July are plotted in green for comparison.

Appendix D: EC flux processing procedure

Two approaches are used to process LES results for EC flux as introduced in Sect. 3.4. The comparison of the two approaches are shown in Fig. D1. The plot shows that most of the times the two approaches generate same results. The small discrepancies in some hours might be due to under-resolved turbulent transport which is missing in the temporal averaging approach.

Figure D1 Comparison of EC flux computed from Eq. (12) and from virtual measurement. Both are using PALM outputs.

Appendix E: Additional information on grid sensitivity study

As stated in the main text, the grid sensitivity study was done on a smaller 2 km by 2 km area in Camden, London. Figure E1 shows the topography of the area in the 10 and 5 m resolution grid. The rasterisation can remove some low buildings and fill in holes in the topography where one single fluid grid is surrounded by solid grids. Figure E2 shows the constant emission patterns applied to the bottom surface of the computational domain. In case of grid in buildings, the emission is applied to the top of the building.

Figure E1 Buildings and terrain heights after rasterisation to 10 and 5 m resolution.

Figure E2 Emission pattern used in the grid sensitivity study. A constant emission grid-map in time is used. The same pattern is used for both 10 and 5 m simulations. In 5 m simulation, the emission grid-map is interpolated bi-linearly from 10 to 5 m resolution.

Appendix F: VPRM parameters

Table F1 lists all the parameters used in the VPRM module in this study. It is a set of parameters fit for Europe (Lian et al., 2021).

Table F1 VPRM parameter values for different vegetation types.

Appendix G: Diurnal cycles of surface CO₂ concentration

Figure G1 shows the diurnal cycles of mean surface CO₂ concentration. This is a complement to the concentration and ABL height analysis in Sect. 4.3.

Figure G1 Diurnal cycles of simulated mean surface level CO₂ concentration for Janurary and July.