Impact of urbanization on gas-phase pollutant concentrations: a regional-scale, model-based analysis of the contributing factors

Urbanization or rural–urban transformation (RUT) represents one of the most important anthropogenic modifications of land use. To account for the impact of such process on air quality, multiple aspects of how this transformation impacts the air have to be accounted for. Here we present a regional-scale numerical model (regional climate models RegCM and WRF coupled to chemistry transport model CAMx) study for present-day conditions (2015–2016) focusing on a range of central European cities and quantify the individual and combined impact of four potential contributors. Apart from the two most studied impacts, i.e., urban emissions and the urban canopy meteorological forcing (UCMF, i.e., the impact of modified meteorological conditions), we also focus on two less studied contributors to the RUT impact on air quality: the impact of modified dry deposition due to transformed land use and the impact of modified biogenic emissions due to urbanization-induced vegetation modifications and changes in meteorological conditions affecting these emissions. To quantify each of these RUT contributors, we performed a cascade of simulations with CAMx driven with both RegCM and WRF wherein each effect was added one by one while we focused on gas-phase key pollutants: nitrogen, sulfur dioxide (NO₂ and SO₂), and ozone (O₃).

The validation of the results using surface observations showed an acceptable match between the modeled and observed annual cycles of monthly pollutant concentrations for NO₂ and O₃, while some discrepancies in the shape of the annual cycle were identified for some of the cities for SO₂, pointing to incorrect representation of the annual emission cycle in the emissions model used. The diurnal cycle of ozone was reasonably captured by the model.

We showed with an ensemble of 19 central European cities that the strongest contributors to the impact of RUT on urban air quality are the urban emissions themselves, resulting in increased concentrations for nitrogen (by 5–7 ppbv on average) and sulfur dioxide (by about 0.5–1 ppbv) as well as decreases for ozone (by about 2 ppbv). The other strongest contributor is the urban canopy meteorological forcing, resulting in decreases in primary pollutants (by about 2 ppbv for NO₂ and 0.2 ppbv for SO₂) and increases in ozone (by about 2 ppbv). Our results showed that they have to be accounted for simultaneously as the impact of urban emissions without considering UCMF can lead to overestimation of the emission impact. Additionally, we quantified two weaker contributors: the effect of modified land use on dry deposition and the effect of modified biogenic emissions. Due to modified dry deposition, summer (winter) NO₂ increases (decreases) by 0.05 (0.02) ppbv, while there is almost no average effect for SO₂ in summer and a 0.04 ppbv decrease in winter is modeled. The impact on ozone is much stronger and reaches a 1.5 ppbv increase on average. Due to modified biogenic emissions, a negligible effect on SO₂ and winter NO₂ is modeled, while for summer NO₂, an increase by about 0.01 ppbv is calculated. For ozone, we found a much larger decreases of 0.5–1 ppbv.

In summary, when analyzing the overall impact of urbanization on air pollution for ozone, the four contributors have the same order of magnitude and none of them should be neglected. For NO₂ and SO₂, the contributions of land-use-induced modifications of dry deposition and modified biogenic emissions have a smaller effect by at least 1 order of magnitude, and the error will thus be small if they are neglected.

1 Introduction

Urbanization represents one of the most important transformations of land use, turning the natural surface into a built-up surface with objects like buildings, streets, and roads. While urban areas only represent less than 1 % a percent of the total Earth surface (Gao and O'Neill, 2020), more than half of the Earth's population already lives in cities (UN, 2018a), and this transformation, which is often called rural–urban transformation (RUT), is still an ongoing process. It is expected that in the upcoming decades, more than 60 % of the population will live in urban areas (UN, 2018b), making research focused on their environmental effects more and more crucial.

It is known that urban areas predominantly affect the atmospheric environment (Folberth et al., 2015), and they act via two primary intrusions that urbanization represents within the natural environment: (i) the introduction of urban land surface replacing rural land surface, causing significant modifications of the meteorological conditions (Oke et al., 2017) and climate (Huszar et al., 2014; Zhao et al., 2017), and (ii) the introduction of a massive emissions source of anthropogenic pollutants perturbing not only local but also regional and global air composition (Lawrence et al., 2007; Timothy and Lawrence, 2009; Im and Kanakidou, 2012; Huszar et al., 2016a).

As for the air quality of urban areas and those surrounding large cities, it is clear that the main driver affecting the concentrations are local urban emissions. Indeed, many studies looked at the perturbation of the atmospheric composition due to solely urban emissions over different scales. For example, Lawrence et al. (2007), Butler and Lawrence (2009), and Stock et al. (2013) investigated the global impact of emissions from megacities, while on regional scales many studies focused on large agglomerations in Europe, like Athens, Istanbul, London, and Paris (e.g., Im et al., 2011a, b; Im and Kanakidou, 2012; Finardi et al., 2014; Skyllakou et al., 2014; Markakis et al., 2015; Hodneborg et al., 2011; Huszar et al., 2016a; Hood et al., 2018), or on large eastern Asian pollution hot spots (Guttikunda et al., 2003, 2005; Tie et al., 2013). These studies show that, not surprisingly, the concentrations of primary pollutants like oxides of nitrogen and sulfur (NO_x, SO₂), volatile organic compounds, and primary aerosols are substantially increased. But on the other hand, urban emissions, due to their high NO_x-to-VOC ratio (VOC: volatile organic compound), can lead to decreases in ozone in the urban cores (e.g., Huszar et al., 2016a). There is further a general consensus in these studies that although air pollution in cities is determined mainly by local sources, a significant fraction of the total concentration is associated with rural sources or sources from other cities (Panagi et al., 2020; Thunis et al., 2021; Huszar et al., 2021).

Urbanization, however, influences the final air pollution in other ways too. One of the most studied aspects of RUT is the modulation of the pollutant concentration due to the meteorological forcing represented by the urban canopy, which includes effects like increased temperatures (urban heat island, UHI) (Oke, 1982; Oke et al., 2017; Karlický et al., 2018, 2020; Sokhi et al., 2022), lower wind speeds (Jacobson et al., 2015; Zha et al., 2019), or elevated boundary layer height along with enhanced vertical eddy diffusion (Ren et al., 2019; Huszar et al., 2020a; M. Wang et al., 2021). Huszar et al. (2020a) introduced the term urban canopy meteorological forcing (UCMF), which represents the forcing that the land surface modified by RUT represents for the physical state of the air above via perturbed exchange of momentum, heat, radiation, and moisture. UCMF is thus a modification of meteorological conditions, which in turn propagates to modifications of pollutant concentrations via modifying the transport, chemical transformation, and deposition of air pollutants. Indeed, Ulpiani (2021) argued that urban pollution has to be studied in connection with UHI and other related meteorological effects. Many other studies looked at the impact of UCMF on air quality and found that the most important parameters in this regard are temperature, turbulence, and wind (Struzewska and Kaminski, 2012; Liao et al., 2014; Kim et al., 2015; Zhu et al., 2017; Zhong et al., 2018; Li et al., 2019; Huszar et al., 2014, 2018a, 2020b), while moisture effects were rather minor (Huszar et al., 2018b). These studies found that these changes led to near-surface decreases in primary pollutant concentrations, while in the case of secondary pollutants (e.g., ozone) increases are encountered either on the surface or at higher levels (Huszar et al., 2018a; Janssen et al., 2017; Yim et al., 2019; Li et al., 2019; Kim et al., 2021; Kang et al., 2022). In other words, besides the urban emission input, UCMF is another factor that contributes to the final urban pollution within the overall process of RUT (Huszar et al., 2021).

Moreover, during urbanization the land use is modified from rural (or natural like forest and grassland) to “urban”, which itself introduces a forcing via a further pathway: contrary to wet deposition, dry deposition velocities (DVs) greatly depend on the land use type, which determines the resistance of the surface and canopy layer (Zhang et al., 2003; Cherin et al., 2015; Hardacre et al., 2021). In urban environments, vegetation is greatly reduced (expressed, for example, in terms of the leaf area index – LAI – reduction). As plants represent a major sink for many gaseous air pollutants (via stomatal uptake), it is clear that over urban areas, this sink is missing or is strongly reduced. For example, based on a later study, over urban land surface the typical DVs of nitrogen dioxide (NO₂) and ozone (O₃) are about half of that above agricultural land like crops (as a typical rural land surface type). Nowak and Dwyer (2007) calculated a net average pollutant removal if additional trees were planted in select US urban areas. Mcdonald-Buller et al. (2001) also showed that ozone and NO₂ are greater in a photochemical model if the land use information supplied contains a higher fraction of urban land use type. In general it seems that, besides other effects, urbanization also leads to increased ozone concentrations due to reduced deposition values (Song et al., 2008; Tao et al., 2015), while dry deposition itself is an important factor determining ozone pollution (Galmarini et al., 2021). However, for some secondary pollutants like HNO₃, H₂SO₄, H₂O₂, HONO, or NH₃, the removal in the case of wet canopies is higher for urban areas than for rural ones (e.g., crops) due to their high solubility and reactivity with solid surfaces (Zhang et al., 2003). Urbanization in the case of these species means higher DVs, leading to a decrease in their concentrations. It is thus clear that the final air pollution caused by urbanization has another contribution represented by the modified (increases and decreases too) dry deposition uptake potential of urban land surface compared to rural and/or natural land surface.

Finally, vegetation not only acts as a sink of pollutants via dry deposition, but it also emits large quantities of biogenic hydrocarbons (biogenic volatile organic compounds – BVOCs; Kesselmeier and Staudt, 1999). Due to their reactivity and potential to form peroxy radicals, they contribute to the formation of tropospheric ozone (Situ et al., 2013; Tagaris et al., 2014). As already mentioned above, during urbanization, the vegetation is strongly reduced, which will result in a decrease in BVOC emissions. Song et al. (2008), for example, showed up to 10 % reductions due to urbanization in Texas. As urban areas are usually VOC-limited environments, reduced BVOC emissions are expected to lead to reduced ozone concentrations (Song et al., 2008). It has to be noted that anthropogenic emissions from urban areas encompass the emissions of VOC compounds, typically of biogenic origin (like isoprene and monoterpenes; Wagner and Kuttler, 2014; Panopoulou et al., 2020). These emissions are probably much smaller than other VOC emissions (Guo et al., 2022) and cannot outweigh the reduction due to reduced vegetation.

Urbanization-induced BVOC emission modifications have a further sub-component acting via modified meteorological conditions in cities. Indeed, urban temperatures are higher than rural ones, and there is an indication that urban cloudiness, at least for European cities, is slightly reduced too (Karlický et al., 2020). These effects have a direct impact on the biochemistry of plants and thus on the quantity of emitted BVOCs as higher temperatures and more solar radiation promote these emissions (Guenther et al., 2006). This means that, due to urbanization, BVOC emissions are suppressed by reducing the vegetation fraction; however, more favorable weather conditions act in an opposite way, making these effects counteracting, although there is an indication that the former (vegetation) effect is dominant (Li et al., 2019).

In summary, urbanization substantially affects air quality, while the final urban pollutant concentration levels are a result of multiple impacts that add to the background (i.e., that without urbanization) air pollution. These include the following.

1. The effect of urban emissions (DEMIS)

2. The effect of the urban canopy meteorological forcing (UCMF) on pollutant transport and chemistry (DMET)

3. The effect of modified dry deposition associated with modified land cover (DLU_D)

4. The effect of modified emissions of biogenic volatile compounds (BVOCs) due to modified land cover (i) and meteorology (ii) (DBVOC)

As seen above, many studies looked at the total impact of urbanization or at some of the individual contributors listed. However, they did not systematically analyze the impact of each one of them. Here we propose a study to uncover the (i) total impact of urbanization (DTOT) and, more importantly, the contribution of (ii) each of the urbanization-related impacts (i.e., DEMIS, DMET, DLU_D, and DBVOC) over a regional-scale domain to present-day urban air pollution levels using coupled regional climate and chemistry transport models applied at moderate 9 km × 9 km horizontal resolution. To achieve this goal, we have to define the reference (or background; not to be confused with “background ozone”, which is a well-defined term) state to which these impacts will be gradually added: rural land use without the effect of UCMF and only rural emissions (urban emissions removed, i.e., those falling within the city administrative boundaries), while present-day land use, emissions, and climate are considered (2015–2016, see below). To evaluate the individual contributors to urbanization as well as their combined effect, we will gradually add each of them to the base simulation in a cascading fashion. To reduce the uncertainty of the results caused by the different geographical and climatic conditions of cities, we perform our analysis for a large ensemble of cities in central Europe: 19 cities in total. Although a similar estimate across several urban areas was made in Huszar et al. (2016a), they focused on the effect of emissions only (which corresponds to DEMIS in our study), while none of the other effects (UCMF; effect of land use on dry deposition and effect of modified biogenic emissions) were considered. Further, it is clear that to some degree the urbanization of one geographic location will affect the air pollution of other urban areas; however, this effect was evaluated to be minor for emissions and UCMF if the cities analyzed are sufficiently far from each other (Huszar et al., 2014, 2016a). Therefore, the selection of cities considers the requirement of sufficient distance from each other.

The study will focus on the key gas-phase pollutants NO₂, O₃, and SO₂. NO₂ is one of the most important primary pollutants in urban environments responsible for reduced air quality and as a precursor for secondary pollutants like ozone or inorganic fine aerosol (Im and Kanakidou, 2012; Stock et al., 2013; Mertens et al., 2020). Ozone is formed in urban plumes when NO_x and VOCs mix together promoted by solar radiation (Xue et al., 2014). Finally, sulfur dioxide is a pollutant originating mainly from fossil fuel combustion in energy production (Guttikunda et al., 2003); although it has undergone significant reduction in European cities during the last decades, it remains of concern, especially in eastern European countries (e.g., in Poland; EEA, 2019).

Of course, urban air pollution is affected not only by local effects. Emissions from other areas (rural or other, even distant cities) constitute a major fraction of urban air pollution (Im and Kanakidou, 2012; Huszar et al., 2016a). Further, background regional air pollution is an important factor that plays a role, e.g., in the urban ozone burden (Yan et al., 2021). Also, UCMF can have regional effects, and the UCMF due to one city can have an effect others (Huszar et al., 2014). However, in this study we are interested in the local effects, i.e., the effect of rural–urban transformation on the local final air pollution, and concerned about the effect of the background atmosphere or the effect from other urban areas.

The study is structured as follows: after the Introduction, the experimental tools (models), their configuration, and the data used are presented. Next, the experiments performed are presented, followed by the Results section. Finally, these aspects are discussed and conclusions are drawn.

2 Methodology

2.1 Models used

The study is based on numerical experiments carried out using regional climate models (RCM) coupled to a chemistry transport model (CTM). To describe the regional climate, two RCMs as meteorological driver are used: RegCM version 4.7 and WRF version 4.0.3. Chemistry was resolved with the chemical transport model CAMx in version 7.10. The decision behind choosing two regional meteorological drivers is to achieve, at least to some degree, more robust results given the fact that the modeled meteorological conditions over cities greatly impacts the chemical concentrations (Ďoubalová et al., 2020; Huszar et al., 2018b).

As the models used and the parameterizations applied are almost identical to those in Huszar et al. (2021), here we list the most important details. RegCM4.7 is a regional-scale climate model with both hydrostatic and non-hydrostatic dynamics (Giorgi et al., 2012). The schemes adopted are the Tiedtke scheme (Tiedtke, 1989) for convection, the Holtslag scheme (HOL; Holtslag et al., 1990) for planetary boundary layer (PBL) parameterization, and the five-class WSM5 moisture scheme (Hong et al., 2004) for microphysics. The atmosphere–biosphere–surface exchange in RegCM was calculated using the Community Land Model (CLM) version 4.5 (Oleson et al., 2013) land surface scheme, and to resolve the urban-scale meteorological phenomena the CLMU module within CLM4.5 is invoked (Oleson et al., 2008, 2010). CLMU considers the traditional canyon geometry approach, meaning that cities are represented as networks of street canyons with specified geometrical and surface parameters (Oke et al., 2017).

WRF (Weather Research and Forecasting Model) is a regional weather prediction and climate model with a detailed description provided by Skamarock et al. (2019). In our modeling setup, the Grell 3D convection scheme (Grell, 1993), the BouLac PBL scheme (Bougeault and Lacarrère, 1989), and the Purdue Lin scheme (Chen and Sun, 2002, PLIN) for microphysics were used. The urban canopy meteorological effects were resolved using the Single-Layer Urban Canopy Model (SLUCM; Kusaka et al., 2001).

For the chemistry simulations the chemistry transport model CAMx version 7.10 (Ramboll, 2020) was used (i.e., we used the most up-to-date version for CAMx available). CAMx is an Eulerian photochemical CTM implementing multiple gas-phase chemistry schemes (Carbon Bond 5 and 6, SAPRC07TC, etc.) with the Carbon Bond 6 revision 5 (CB6r5) scheme used in this study. CB6r5 includes updates to chemical reaction data from IUPAC (IUPAC, 2019) and NASA (Burkholder et al., 2019) for inorganic and simple organic species important for the formation of ozone. Apart from the inclusion of the CB6r5 mechanism, this version of CAMx includes important modifications of secondary aerosol formation via oxidation of VOCs, which can have feedbacks on the total VOC and thus ozone concentrations. To complete the atmospheric chemistry with aerosol physics, a static two-mode approach was considered. The ISORROPIA thermodynamic equilibrium model (Nenes et al., 1998) was invoked for the secondary inorganic aerosol formation. Secondary organic aerosol (SOA) was partitioned from its gas-phase precursors using the SOAP equilibrium scheme (Strader et al., 1999). For wet deposition, the Seinfeld and Pandis (1998) scheme is used, while dry deposition is treated using the Zhang et al. (2003) method. The Zhang method incorporates a three-resistance equation for deposition velocity (DV) incorporating the aerodynamic resistance (r_a), the above-canopy quasi-laminar sublayer resistance (r_b), and the overall canopy resistance (r_c), while DV is calculated as the reciprocal of the sum of these resistances. An important component of r_c is the resistance represented by vegetation with the so called in-canopy aerodynamic, stomatal, mesophyll, and cuticle resistances. Over urban surfaces or any other non-vegetated surfaces, these are not defined; however, to use the same equations for such land use categories in the dry deposition model, very large values are applied (e.g., 10²⁵ sm⁻¹). Regarding the aerodynamic resistance representing the bulk transport trough the lowest model layer via turbulent diffusion, its magnitude depends on the intensity of turbulence, which in turn depends on wind speed, surface roughness, near-surface temperature lapse rate, and solar insolation. Over urban areas these are also strongly modulated, implying a strong influence on the deposition velocities. Finally, the quasi-laminar sublayer (or boundary) resistance r_b represents molecular diffusion through the thin layer of air directly in contact with the surface, and this is assumed to be a function of the molecular diffusivity of each pollutant regardless of the surface it is deposited on. In this dry deposition model, the aerodynamic resistance has a strong dependence on the temperature via decreased stability near the surface and thus more efficient turbulent diffusion towards the surface (Louis, 1979). Also, the stomatal resistance decreases with higher temperatures due to wider stomata (Zhang et al., 2002); this holds up to a threshold maximum temperature at which stomata suddenly close. These temperatures are usually, however, not reached in the climate of the region. Therefore, increased dry deposition velocities are expected as the result of increased temperatures in urban areas.

A meteorological preprocessor is used to convert the RegCM and WRF meteorological data into model-ready driving data for CAMx: for the WRF, it was the wrfcamx preprocessor, which is provided along with the CAMx code at https://www.camx.com/download/support-software/ (last access: 27 September 2022), while for RegCM, the RegCM2CAMx interface was applied (Huszar et al., 2012). The vertical eddy diffusion coefficients (K_v) are diagnosed from the available meteorological data on RegCM and WRF output using the CMAQ diagnostic approach (Byun and Ching, 1999). Temperature, pressure, humidity, and cloud–rain water content are defined at cell centers along with pollutant concentration, and CAMx considers them to be grid cell average conditions. On the other hand, wind and diffusion coefficients are carried at cell interfaces to describe the mass transfer across each cell face. Coupling between CAMx and the driving models is offline, which implies that no feedbacks of the pollutant concentrations on WRF/RegCM radiation and microphysical processes were considered. Indeed, Huszar et al. (2016b) showed with 10-year-long simulations that the long-term radiative effects of urban pollutant emissions and secondarily formed pollutants (like ozone) are rather small, which justified this choice.

2.2 Model setup and data

Model simulations were performed over identical domains (parent and nested ones) and for an identical period as in Huszar et al. (2021), i.e., the years 2015–2016 with 9, 3, and 1 km horizontal resolution centered over the Czech capital, Prague (50.075^∘ N, 14.44^∘ E; Lambert conic conformal projection). In the vertical, the model grid has 40 layers in both meteorological driving models. The thickness of the lowermost layer is about 30 m, and the top of the model's atmosphere reaches 5 hPa (about 36 km). The simulated time period is December 2014–December 2016 (the first month is used as spin-up). Tie et al. (2010) argued that the ratio of the diameter of the analyzed city to model resolution should be at least 6 : 1, which means that in our case, a 6 km or smaller horizontal grid step should be used to resolve the impact of urbanization for the cities chosen (see below). For Prague, which is modeled at 1 km, this is fulfilled. Other cities outside the inner 1 km nested domain are treated at coarser resolution, but many studies found that the impact of the resolution of emissions and models on urban species concentrations is rather small: e.g., Hodneborg et al. (2011) showed that coarse resolution can lead to higher ozone modeled by around 10 %, while Markakis et al. (2015) found only moderate sensitivity (8 %) to model resolution. Similarly, Y. Wang et al. (2021) showed that ozone production is reduced when high resolution is applied, but the reduction is only about 8 % for ozone.

The ERA-Interim reanalysis (Simmons et al., 2010) is used as forcing data. The 3 and 1 km domains are then driven by the corresponding parent domains with one-way nesting. Chemical boundary conditions are based on the CAM-chem global model data (Buchholz et al., 2019; Emmons et al., 2020). Land use data (for both climate models and for the dry deposition scheme in CAMx) were derived from the high-resolution (100 m) CORINE CLC 2012 land cover data (https://land.copernicus.eu/pan-european/corine-land-cover, last access: 8 August 2022) as well as from the United States Geological Survey (USGS) database for grid cells with no information from CORINE. In RegCM, fractional land use is considered, while in WRF, each grid cell is attributed the dominant land use, which brings some accounting for the uncertainty related to the urban land cover representation. This means that due to the fact that the land use is represented differently in WRF and RegCM, partly urbanized surfaces and their effects can be differently accounted for in the two models.

The European CAMS (Copernicus Atmosphere Monitoring Service) version CAMS-REG-APv1.1 inventory (Regional Atmospheric Pollutants; Granier et al., 2019) for the year 2015 was used as anthropogenic emission data for areas outside the Czech Republic. There, high-resolution national data were adopted: the Register of Emissions and Air Pollution Sources (REZZO) dataset issued by the Czech Hydrometeorological Institute (https://www.chmi.cz, last access: 27 September 2022) and the ATEM Traffic Emissions dataset provided by ATEM (Ateliér ekologických modelů – Studio of Ecological Models; https://www.atem.cz, last access: 27 September 2022) were used. These data provide activity-based (SNAP – Selected Nomenclature for sources of Air Pollution) annual emission totals of oxides of nitrogen (NO_x), volatile organic compounds (VOCs), sulfur dioxide (SO₂), carbon monoxide (CO), PM_2.5 and PM₁₀ (particles with diameter less than 2.5 and 10 µm), and ammonia (NH₃). CAMS data are defined on a regular Cartesian lat–long grid, while the Czech datasets are provided as area, line (for road transportation), or point sources (in the case of area sources these are usually irregular shapes corresponding to counties with resolution from a few tens of meters to 1–2 km). The Flexible Universal Processor for Modeling Emissions (FUME) emission model (http://fume-ep.org/, last access: 27 September 2022; Benešová et al., 2018) is used to preprocess the mentioned emission inventories to CTM-ready emission files, including preprocessing the raw input files, the spatial remapping of the data into the model grid, chemical speciation, and time disaggregation from annual to hourly emissions. Speciation factors and time disaggregation profiles were taken from Passant (2002) and van der Gon et al. (2011), respectively. The temporal factors contain activity-sector-specific monthly, weekly, and hourly factors used to decompose the annual totals into hourly emissions. Geographic dependence is not considered here; however, the time zone information and the associated time shifts are accounted for.

Emissions of biogenic origin are calculated offline using MEGANv2.1 (Model of Emissions of Gases and Aerosols from Nature version 2.1) (Guenther et al., 2012) based on RegCM and WRF meteorology (temperature, shortwave radiation, humidity, soil moisture). The necessary input for MEGAN including leaf area index data (its annual cycle), plant functional types, and emission potentials of different plant types are not part of the CORINE land use data and were derived independently from Sindelarova et al. (2014, 2022). It has to be mentioned here that along with the calculation of biogenic VOC data, MEGAN also calculates the fluxes of soil biogenic NO (nitrogen monoxide) emissions as a result of bacterial activity in soil according to Yienger and Levy (1995). As these emissions are a function of LAI and meteorological conditions, part of the DBVOC impact will be composed of soil NO_x emissions modifications. Although not presented here, in our experiments the soil NO_x emissions are about 2 orders of magnitude smaller compared to the BVOC emissions, and their effect is expected to be much smaller including the effect of their urbanization-induced modifications. It also has to be stressed that BVOC emissions are strongly temperature-dependent, while higher temperatures trigger stronger emissions. In this regard, urbanization-induced temperature enhancement is expected to lead to stronger BVOC fluxes. Wildfire emissions can potentially be episodically important and can significantly contribute to levels of gaseous pollutants like NO_x and CO as well as to improving overall model performance (Lazaridis et al., 2008); they are, however, significant mainly over southern Europe and the Mediterranean and not over our focus area (central Europe). Moreover, wildfire emissions normally do not occur in urban areas and therefore do not contribute to the impact of urbanization.

A key task was to isolate the emissions originating from urban areas (see further for details about the chosen cities). In this regard, urban areas were identified based on the administrative boundaries of chosen cities. We used the GADM public database (https://gadm.org, last access: 27 September 2022) for their definition. While masking of inventory emissions based on the GADM shapes corresponding to cities, it had to be ensured that the partition between the “city” and “non-city” emission segments (inventory grid cells or irregular shapes in the case of Czech emissions) over the city boundary were correctly calculated. For this purpose, the masking capability of FUME was adopted.

The cities chosen in the analysis are Berlin, Brussels, Budapest, Cluj-Napoca, Cologne, Frankfurt, Hamburg, Krakow, Lodz, Lyon, Milan, Munich, Prague, Torino, Vienna, Warsaw, Wroclaw, Zagreb, and Zurich. They are also highlighted in Fig. 1 including the 9 km domain terrain elevation. The choice of the cities used the same criteria as in Huszar et al. (2021): the size of the city comparable to one 9 km × 9 km grid cell, sufficient distance between cities to eliminate inter-city influences, minimal orographic variability to reduce orographic effects (Ganbat et al., 2015), and no coastal cities to eliminate the effect of asymmetric land use like the sea breeze effect (Ribeiro et al., 2018). Although strict emission control policies, these cities are still often burdened with high air pollution for pollutants as NO₂ and O₃ (EEA, 2019; Khomenko et al., 2021; Sokhi et al., 2022).

Figure 1 The 9 km × 9 km resolution model domain and the resolved terrain in meters including the cities analyzed in the study (red squares).

2.3 Model simulations

The study intends to evaluate the urbanization impact on air quality, while we attempted to decompose the total impact into individual contributors listed in the Introduction. This requires performing a series of model experiments with individual effects added gradually one by one to a reference state to end up with the real situation corresponding to full urbanization. In Huszar et al. (2018a, b) we performed a similar decomposition for the urban-induced meteorological effects (i.e., the UCMF) and their impact on air quality. Here we adopt this approach, but it will not concern the UCMF solely as in the mentioned works but the entire impact of urbanization, while UCMF will be treated as one effect.

The simulations performed are summarized in Tables 1 and 2 for RCM and the underlying CTM simulations, respectively. A pair of simulations was performed with both RegCM and WRF with (“Urban”) and without (“Nourban”) considering urban land surface. In the latter case, land use was replaced by “crops” as the most common rural land use type in the region analyzed. While all RegCM simulations and CAMx simulations driven by RegCM were performed on nested domains (9, 3, and 1 km), the WRF simulations and CAMx ones driven by WRF were done only over the parent 9 km domain as WRF served as a complementary model to account for the uncertainty in the driving meteorology, especially with regard to UCMF; note that the urban canopy model is different in WRF than in RegCM.

Table 1 The list of RCM simulations performed.

^a Information on whether urban land surface was considered. “Nourban” means replacing the urban land surface by crops. ^b Simulation performed in a nested way at 9, 3, and 1 km horizontal grid resolution.

Table 2 The list of CTM simulations performed with information on the effects considered. The “driving meteorology” and “driving meteorology (BVOC)” columns correspond to the information from Table 1 above, i.e., which RCM simulation from the “Urban”–“Nourban” pair was used.

^∗ Information on whether the meteorology driving MEGAN accounted for the UCMF.

As for the CAMx simulations, they differ based on the inclusion of urbanized–rural land surface, the UCMF (acting on both atmospheric chemistry in general and on BVOC fluxes), and the urban emissions. In this regard, we performed six experiments summarized in Table 2. The reference experiment called ENNrrN represents the hypothetical background state without urban emissions and with the urban land surface replaced by rural land surface in RCMs and the CTM as well as in the BVOC model (MEGAN). The reference simulation is not to be confused with the preindustrial state: we assumed current climate (GHG – greenhouse gas – concentration) and current background large-scale chemical concentrations. In the next experiment, ENYrrN, only the urban emissions are considered (turned on). In the third experiment, ENYurN, the urban land use was turned on for the dry deposition in CAMx. In the fourth experiment, ENYuuN, the urban land use is also turned on for the biogenic emission model. In the fifth experiment, ENYuuU, both the urban land use and the UCMF are accounted for in the biogenic emissions model, and finally, in the sixth experiment, EUYuuU, all the urbanization-related effects are considered, representing the most realistic case.

In the first experiment wherein urban emissions are disregarded, we removed urban emissions only for the 19 cities chosen for the analysis. For the effect of rural–urban land use transformation on meteorological conditions, dry deposition, and biogenic emissions, we replaced the urban land by rural land over the entire domain (i.e., not only for the cities chosen). It is clear that this has an effect on the background level of air pollutants not only at local urban levels, but the effect is probably much smaller than local effects as (1) emissions from these areas were still considered, and (2) the urban meteorological effects from these (minor) urban areas have a rather small influence on air pollutants as the UCMF is also small (see, e.g., Huszar et al., 2014). In Fig. 1, we plot the model orography and the analyzed cities as red squares. For urban land use information used in our study, please refer to Karlický et al. (2020) (Figs. 1 and 2), who used identical land use as in our study.

Mathematically, with respect to the rural–urban transformation (RUT), the concentration c_i of a pollutant i for a chosen city is given by

$$\begin{array}{}\text{(1)}& c_i=c_{i,\mathrm{rural}}+\mathrm{\Delta }{c}_{i,\mathrm{RUT}},\end{array}$$

where c_i,rural is the average concentration before RUT and Δc_i,RUT is the total impact of urbanization.

In this study, we are concerned about the contributors to Δc_i,RUT (regardless of their sign):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{c}_{i,\mathrm{RUT}}=\mathrm{\Delta }{c}_{i,\mathrm{EMIS}}+\mathrm{\Delta }{c}_{i,\mathrm{MET}}+\mathrm{\Delta }{c}_{i,\mathrm{LU}\mathrm{\_}\mathrm{D}}+\mathrm{\Delta }{c}_{i,\mathrm{BVOC}},\end{array}$$

where Δc_i,EMIS, Δc_i,MET, Δc_{i,LU_D}, and Δc_i,BVOC are the impacts of urban emissions, the impact of the urban canopy meteorological forcing, the impact of modified land use on dry deposition, and the impact of modifications of BVOC emissions, denoted above as DEMIS, DMET, DLU_D, and DBVOC. The Δc_i,BVOC impact can further be decomposed into the part caused by modified land cover (reduced vegetation in terms of changes in leaf area index – LAI, DBVOC_L) and modified meteorological conditions (DBVOC_M):

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{c}_{i,\mathrm{BVOC}}=\mathrm{\Delta }{c}_{i,\mathrm{BVOC}\mathrm{\_}\mathrm{L}}+\mathrm{\Delta }{c}_{i,\mathrm{BVOC}\mathrm{\_}\mathrm{M}}.\end{array}$$

These impacts will be calculated from the experiments listed in Table 2 in the following way (the experiment number is shown in parentheses).

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{\Delta }{c}_{i,\mathrm{RUT}}=\mathrm{EUYuuU}\left(\mathrm{6}\right)-\mathrm{ENNrrN}\left(\mathrm{1}\right)\text{(5)}& {\displaystyle }\mathrm{\Delta }{c}_{i,\mathrm{EMIS}}=\mathrm{ENYrrN}\left(\mathrm{2}\right)-\mathrm{ENNrrN}\left(\mathrm{1}\right)\end{array}$$

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathrm{\Delta }{c}_{i,\mathrm{MET}}=\mathrm{EUYuuU}\left(\mathrm{6}\right)-\mathrm{ENYuuU}\left(\mathrm{5}\right)\text{(7)}& {\displaystyle }\mathrm{\Delta }{c}_{i,\mathrm{LU}\mathrm{\_}\mathrm{D}}=\mathrm{ENYurN}\left(\mathrm{3}\right)-\mathrm{ENYrrN}\left(\mathrm{2}\right)\text{(8)}& {\displaystyle }\mathrm{\Delta }{c}_{i,\mathrm{BVOC}}=\mathrm{ENYuuU}\left(\mathrm{5}\right)-\mathrm{ENYurN}\left(\mathrm{3}\right)\end{array}$$

It has to be noted that in reality these effects act simultaneously and feedbacks are present between them, so their effects are not additive. The way we calculated the individual impacts (contributors), however, allows us to consider them to be additive; i.e., their sum is the total impact of urbanization.

3 Results

3.1 Validation

This model configuration (same input data, same domain) underwent a detailed validation including both meteorology and air quality in Huszar et al. (2020b, 2021). However, due to the fact that in our CAMx simulations, the newer version 7.10 was used instead of version 6.50, and instead of the CB5 we used the newer CB6 chemistry mechanism, we provide a brief account for validation. For comparison with observations, AirBase European air quality data (https://discomap.eea.europa.eu/map/fme/AirQualityExport.htm, last access: 27 September 2022) for the modeled years were used, while all urban- and suburban-type background stations (AirQualityStationArea indicates urban and/or suburban, and AirQualityStationType indicates background) were used from a subset of the analyzed cities. This sub-selection considered the largest cities from the total of 19 where sufficient numbers of stations were available.

Figure 2 Comparison of modeled (blue – RegCM/CAMx; green – WRF/CAMx) and measured (red; AirBase data) urban and suburban average monthly concentrations of NO₂ (first column) and O₃ (second column), the average JJA diurnal cycle of O₃ (third column), and the average monthly concentrations of SO₂ (fourth column) for 8 different cities selected from the total 19 considered in the study, namely Berlin, Budapest, Milan, Munich, Prague, Zurich, Vienna, and Warsaw (units: µg m⁻³). Data are averaged across all available urban- and suburban-type background stations within the chosen city. There are no data for SO₂ in Munich as no corresponding measuring station was available.

In Fig. 2, the comparison of average monthly means of modeled and measured concentrations of the three analyzed pollutants is shown, while the full experiment (EUYuuU) was used, which represents the real case. For ozone, we also included the average summer diurnal cycle, as daily peak values are more important for this pollutant than the averages values. For Prague, results are taken from the 1 km nested domain; otherwise, they are extracted from the 9 km regional domain. For NO₂, there is a generally acceptable match between the model and observations, with model biases up to 10 µg m⁻³. While concentrations from January to April are usually underestimated, during summer, CAMx generates a positive bias, except in Berlin, where there is an underestimation of NO₂ from March to September and an overestimation during the rest of the year. During late autumn the model bias is usually negative, with large differences between cities. The WRF-driven CAMx concentrations are usually lower then for RegCM/CAMx during summer, which means higher model bias (up to −10 µg m⁻³). During winter, WRF/CAMx gives higher concentrations than RegCM/CAMx, resulting in a smaller bias.

For O₃, both RegCM/CAMx and WRF/CAMx capture the annual cycle well, with some overestimation of concentrations during late spring (by about 10–20 µg m⁻³) and an underestimation during late summer (by a similar magnitude). In general, these two simulations are very similar. During winter, there is a small negative (up to −10 µg m⁻³) model bias present. The summer diurnal cycles show good agreement in the basic measured pattern, including the timing of the maximum. The maximum values are sometimes underestimated (by up to 20 µg m⁻³ for some cities), especially for the RegCM-driven runs. In general, WRF meteorology causes higher simulated maximum ozone. During nighttime, ozone is often overestimated by around 10–20 µg m⁻³, and clearly, the WRF-driven run performs better during this part of the day.

Regarding SO₂, the model fails to capture the annual cycle well. During winter, both models usually underestimate the concentrations by up to 1–2 µg m⁻³ except in Budapest and Berlin, where overestimation occurs in the model. During summer, measured SO₂ concentrations are usually smaller and the models somewhat reflect this fact, but large biases are still present and the models are unable to correctly capture the annual cycle for some cities (e.g., Budapest and Vienna).

3.2 The overall impact of individual contributors to RUT

Firstly, we evaluated the impact of individual contributors to the RUT as well as the total impact in terms of 2015–2016 DJF and JJA averages (in the case of ozone as a summer average only), averaged across the chosen cities and the model ensemble (i.e., average of the RegCM- and WRF-driven runs). Values are taken from the grid box covering the center of a particular city. The results are shown in Fig. 3 as box plots showing the first and third quartiles as well as the median values and the minimum and maximum. The analysis showed (as expected) that from the four contributors to RUT, two are much stronger than the other two. Therefore, in the plots, we separated them from the minors ones (including the total impact).

Figure 3 The 2015–2016 DJF and JJA averaged impact of each contributor to the rural–urban transformation including the total impact averaged over all chosen cities for NO₂, SO₂, and O₃ in parts per billion by volume (ppbv). In the case of ozone, only the summer averaged MDA8 (maximum daily 8 h average) is shown. The box plots show the 25 % to 75 % quantiles including the minimum and maximum value across all cities. The red line shows the median. Values are taken from the model grid cell that covers the city center. Panels (a)–(c) show the two main contributors including the total impact (DEMIS, DMET, and DTOT), while panels (d)–(f) show the minor contributors (DLU_D and DBVOC).

For all three gas-phase pollutants, the impact of emissions (DEMIS) is the largest in magnitude in both seasons. For NO₂ it ranges (i.e., the 25th to 75th percentile) from 4 to about 8 ppbv and from 5 to 10 ppbv for JJA and DJF, respectively. For SO₂, the numbers are somewhat smaller, as cities, at least in the region in focus, are not such strong SO₂ emitters (compared to NO₂): an increase by 0.4 to 1.5 and 0.8 to 1.6 ppbv for JJA and DJF, respectively, is seen. For O₃, the impact on the summer maximum daily 8 h average concentration (MDA8) is characterized by a decrease due to titration (as expected) by 3 to 6 ppbv.

The impact of the urban canopy meteorological forcing (DMET) is characterized by a decrease for the primary pollutants: for NO₂, the decrease is usually between 1 and 6 ppbv for JJA and between 1 and 3 ppbv in DJF with the maximum surpassing zero, meaning that in some cities, a slight increase was modeled. In the case of SO₂, the impact of UCMF is smaller, with up to a 0.6 and 0.4 ppbv decrease in JJA and DJF, respectively. For O₃, the impact is an increase of about 2 ppbv.

In the case of minor contributors, the impact of BVOCs is considerable for ozone only, and, as expected, for the other pollutants it acts as a minor modulator of the overall chemistry (e.g., influencing the hydroxyl budget), therefore having a very small impact. In the case of NO₂, the impact is a slight increase by around 0.01 ppbv in JJA and negligible in winter. For SO₂, it is near zero in both seasons. For ozone, which is directly influenced by biogenic emissions, the impact is a decrease by around 0.4 to 1 ppbv as a JJA average. The impact of modified dry deposition due to urbanized land use is characterized by an increase (0.02 to 0.08 ppbv) for NO₂ in JJA and an opposite impact in winter (around 0.01 to 0.04 ppbv decrease). For SO₂, the impact in summer can be both negative and positive (from −0.01 to 0.01 ppbv) with a near-zero average. In winter, there is a decrease by about 0.02 to 0.07 ppbv. For ozone, the impact of land use change is an increase between 1 and 2 ppbv.

Finally, the total impact is an increase for all pollutants and quantities: for NO₂ it is about 1–5 ppbv in JJA and 4–8 ppbv in DJF, and for SO₂ it ranges from 0 to 1 ppbv in JJA and about 0.5 to 1 ppbv in DJF. For JJA MDA8 ozone, the total impact is characterized by an increase up to 2 ppbv.

3.3 The spatial distribution of the impacts

The box plots presented above give an overview of the averaged impact across all the cities including the distribution around the median value. To obtain spatially resolved information, we also plotted the 2D distribution of the individual contributors here.

3.3.1 The impact of urban emissions (DEMIS)

In Fig. 4 the DJF and JJA average spatial impact of urban emissions (DEMIS) on the near-surface concentrations of NO₂, SO₂, and O₃ is shown.

Figure 4 The spatial distribution of the 2015–2016 average emission impact (DEMIS) for NO₂ DJF and JJA (first and second row), SO₂ DJF and JJA (third and fourth row), and JJA MDA8 O₃ (fifth row). Columns represent the results from the 1 km RegCM/CAMx (detail of Prague), 9 km RegCM/CAMx, and 9 km WRF/CAMx simulations. Units are in parts per billion by volume (ppbv).

In the case of NO₂, the impact reaches 4–6 ppbv in the core of the cities and remains high over surrounding areas (up to 0.5 ppbv over large areas in DJF, especially in the WRF-driven simulations). In summer, the spatial extent of the emission impact is smaller: below 0.1 ppbv over rural areas. The result from Prague at high resolution reveals that the high emission impact is concentrated in the very center of the city (reaching 4–6 ppbv).

For SO₂, there is a larger spread between cities with large contributions over Poland reaching 6 ppbv (in both seasons), while for other cities, the contribution is smaller: up to 2–3 ppbv. The contribution over rural areas is large in Poland (up to 0.2 ppbv) but remains below 0.1 ppbv in other regions. The impact of emissions over Prague reaches 1 ppbv in DJF with contributions up to 0.1 ppbv in its vicinity. In summer, due to low emissions (SO₂ is emitted largely by heating) the contributions are very small, reaching 0.5 ppbv in some hot spots within the city.

Ozone is usually titrated in city centers, which corresponds to the impact of emissions on its concentrations. They decreased over cities by up to 3–4 ppbv, while further from cities, where urban NO_x mixes with rural emissions, an ozone increase occurs of up to 1 ppbv as MDA8. Over Prague, the decrease is limited to the city area. Over its vicinity, the impact becomes positive with a 0.5–1 ppbv increase (similarly as seen for other cities).

3.3.2 The impact of modified meteorological conditions (DMET)

Figure 5 presents the DJF and JJA average spatial impact of the urban canopy meteorological forcing (DMET) on the near-surface concentrations of NO₂, SO₂, and O₃.

Figure 5 The spatial distribution of the 2015–2016 average impact of the urban canopy meteorological forcing (UCMF) DMET on NO₂ DJF and JJA (first and second row), SO₂ DJF and JJA (third and fourth row), and JJA MDA8 O₃ (fifth row). Columns represent the results from the 1 km RegCM/CAMx (detail of Prague), 9 km RegCM/CAMx, and 9 km WRF/CAMx simulations. Units are in parts per billion by volume (ppbv).

In the case of NO₂, while for RegCM/CAMx the impact is usually characterized by a decrease by 1–2 ppbv with some urban areas even showing an increase (up to 2 ppbv), for WRF/CAMx a clear decrease occurs of up to 3 ppbv. For Prague, the highest decreases are modeled in the city center, reaching about 3 ppbv in summer and about 2 ppbv in winter. In general, the winter impact is comparable to summer (slightly stronger for WRF/CAMx).

For SO₂ the impact is weaker and constitutes both decreases (in cities) and increases (over their vicinity) with changes in the range −3 to 3 ppbv. In the case of WRF/CAMx simulations, the impact is more straightforward with the decrease dominating, reaching 3 ppbv in both seasons.

Finally, summer MDA8 ozone increases due to UCMF by up to 2–3 ppbv over cities, while over rural areas, a slight decrease is modeled of up to 1 ppbv appearing in the RegCM/CAMx simulation. Over Prague, the largest increases are modeled in the city center, reaching 2–3 ppbv.

3.3.3 The impact of dry deposition modifications (DLUC_D)

The impact of urban land cover via modified dry deposition (DLU_D) is plotted in Fig. 6. In general, the impacts are much smaller for NO₂ and SO₂ than seen for the emission or the UCMF impact for these pollutants above. For NO₂, the DJF and JJA impacts differ in sign (in accordance with the box plots seen in Fig. 3) and the spatial distribution is somewhat different in WRF/CAMx than in RegCM/CAMx. In DJF, NO₂ concentrations decreased over cities by up to 0.04 ppbv, with some higher decreases over Italy (Milan) of up to 0.1 ppbv. In the WRF-driven experiment, some increases over the Benelux states are also seen, reaching 0.06 ppbv. For Prague, the decrease is maximal in the city center reaching 0.04 ppbv. During JJA, the DLU_D impact is positive, reaching 0.1 ppbv in both models with some slight decreases around Milan. Over Prague, the increase is even stronger and exceeds 0.05 ppbv.

Figure 6 The spatial distribution of the 2015–2016 average impact of the urban land cover via dry deposition modifications (DLU_D) on NO₂ DJF and JJA (first and second row), SO₂ DJF and JJA (third and fourth row), and JJA MDA8 O₃ (fifth row). Columns represent the results from the 1 km RegCM/CAMx (detail of Prague), 9 km RegCM/CAMx, and 9 km WRF/CAMx simulations. Units are in parts per billion by volume (ppbv).

For SO₂, there are clear decreases modeled during DJF, reaching 0.1–0.2 ppbv over city centers. The impacts are slightly stronger in WRF-driven CAMx runs and are about −0.03 ppbv over Prague's center. During JJA, the SO₂ response is very small and positive in the RegCM/CAMx experiments with up to a 0.1 ppbv increase in some city centers, especially over eastern Europe where SO₂ emissions are higher. Decreases similar to the DJF impact remained in the WRF/CAMx simulation. Over Prague, almost zero impact is modeled (between −0.01 and 0.01 ppbv with some positive impact around strong point sources north from the city).

A much stronger response to changes in dry deposition is modeled for summer O₃, with a clear increase reaching 2 ppbv in city centers and being high over rural areas too (up to 1 ppbv increase). Over Prague, the increase is usually 1.5–2 ppbv, exceeding 2 ppbv in the very core of the city.

Figure 7 The spatial distribution of the 2015–2016 average impact of the urban land cover on deposition velocities of NO₂ (a, b), SO₂ (c, d), and O₃ (e, f) for DJF (a, c, e) and JJA (b, d, f) (mm s⁻¹) for the RegCM-driven 9 km CAMx simulations.

To facilitate the interpretation of the simulated responses of concentrations to DLU_D, we also mapped the geographical distribution of the DLU_D impact on the deposition velocities (DV is standardly provided in CAMx output), as seen in Fig. 7 taken from the RegCM-driven CAMx simulations (and not showing the Prague 1 km case). For NO₂, dry deposition velocities decreased by around 0.2 mm s⁻¹ in DJF, and a stronger decrease, reaching −0.6 mm s⁻¹, in city centers is modeled in JJA. For SO₂, the DJF and JJA maps differ in sign. For winter, deposition velocities increased in cities by up to 0.4–0.6 mm s⁻¹, while during summer, similar decreases are simulated compared to NO₂ (around −0.4 to −0.6 mm s⁻¹). For O₃, both seasons are characterized by decreases: by around 0.2 mm s⁻¹ in DJF with stronger decreases in city centers in JJA, reaching −1.5 mm s⁻¹. The WRF/CAMx impacts are very similar and are not shown here.

3.3.4 The impact of biogenic emissions (DBVOC)

The urbanization-induced changes in BVOC emissions (via reduced vegetation cover and modified temperatures; DBVOC) and their consequent effect on summer ozone and NO₂ concentrations are plotted in Fig. 8. As BVOC emissions are of minor importance in winter and the effect on SO₂ is almost zero, we show only the summer impacts for these two pollutants. For NO₂, the DBVOC impact results in increases usually up to 0.06 ppbv, while much stronger increases are modeled over northern Italy (around Milan) of around 0.1 ppbv in both RegCM- and WRF-driven simulations. For Prague, the maximum increase is 0.2–0.3 ppbv.

Figure 8 The spatial distribution of the 2015–2016 average impact of urbanization due to modifications of biogenic emissions (DBVOC) on JJA average NO₂ (a–c) and summer MDA8 O₃ (d–f). Columns represent the results from the 1 km RegCM/CAMx (detail of Prague), 9 km RegCM/CAMx, and 9 km WRF/CAMx simulations. Units are in parts per billion by volume (ppbv).

For O₃, decreases are modeled reaching −1 ppbv over many cities and reaching −0.2 to −0.5 ppbv over rural areas (in the RegCM-driven experiment). A stronger decrease is modeled (again) over northern Italy, with up to −2 ppbv over Milan. Prague is characterized by a decrease usually between −0.5 and −1 ppbv.

Figure 9 The absolute (upper row; units mol km⁻² h⁻¹) and relative change (lower row; units in %) in 2015–2016 JJA averaged isoprene (ISOP) emissions decomposed into the part caused by reduced vegetation (via leaf area index; DBVOC_L) and the part caused by modified meteorology (DBVOC_M). The first and third columns show the change due to DBVOC_L taking the rural (NU) and urban (U) meteorological conditions as a reference, respectively. In the second and fourth columns, the changes due to urban meteorological effects (UCMF) are shown (DBVOC_M) taking the rural (NU) and urban (U) vegetation as a reference, respectively.

The impacts presented above are the result of modified BVOC emissions; therefore, we also plotted the summer changes in isoprene (ISOP) as a major component of such emissions. As changes in these emissions are the result of two components constituted of vegetation change via LAI change (DBVOC_L) and modification of meteorological conditions (the UCMF; denoted DBVOC_M in the Introduction), we plotted the two contributors separately in Fig. 9 as absolute and relative change. We were also interested in whether the reference with respect to which the change is calculated matters. In others words, what is the difference between the DBVOC_L calculated with rural (NOURBAN, see Table 1) meteorology and DBVOC_L calculated with urban meteorology (URBAN in Table 1)? Similarly for DBVOC_M, it is calculated with both rural LAI and that adapted for urban conditions. The impact of vegetation change is an expected decrease in isoprene emissions by up to 15 mol km⁻² h⁻¹, with higher values over the southern part of the domain, often representing a 80 %–90 % decrease in relative numbers, especially for larger and dense urban areas like Milan (Italy). For smaller urban areas the decrease is around −5 % to −20 % (many of the grid cells are only partly covered by urban areas, so the emission decrease is correspondingly small). As seen from the figure, the changes calculated at rural and urban meteorological conditions are very similar (the case with urban meteorology is slightly higher). Regarding the isoprene emission modifications due to UCMF, they are usually much smaller (usually less than 0.05 mol km⁻² h⁻¹ or less than 0.5 % in relative numbers). At some urban areas over Germany as well as over northern Italy and southern France, the change can reach 0.4 to 0.6 mol km⁻² h⁻¹, peaking at 1–2 mol km⁻² h⁻¹ over Italian urban areas, representing a 5 %–10 % relative increase. DBVOC_M is somewhat smaller if calculated with urban land cover, which is expected as the strongest meteorological modifications due to UCMF are over cities, but in this case they affect a non-vegetated surface, which implies smaller effects. In summary, the BVOC emission changes associated with vegetation change are much more important than the modifications due to UCMF.

3.4 The diurnal variation of the impacts

Urban emissions have a strong diurnal cycle caused by the typical cycle of human activities during the day. Moreover, the urban-land-surface-triggered meteorological modifications (UCMF) also have a strong diurnal pattern; e.g., temperature is impacted most during night, and the wind impacts and turbulence modifications are the strongest during noon. (Huszar et al., 2018a, 2020a). Thus, it is clear that the individual contributors to RUT analyzed here are expected to also have a diurnal cycle. Figure 10 presents these cycles for the four contributors and three analyzed pollutants.

Figure 10 The 2015–2016 average diurnal cycle of the individual contributors to RUT for NO₂ (a, d), SO₂ (b, e), and O₃ (c, f) as a DJF (a–c) and JJA (d–f) average. The brown and blue lines stand for the two stronger contributors (DEMIS and DMET, left y axis), while red and green stand for the weaker contributors (DLU_D and DBVOC, right y axis). Units are in parts per billion by volume (ppbv). Times are in UTC (i.e., the local time is +2 h in JJA and +1 h in DJF).

For NO₂ the diurnal pattern for the emissions impact (DEMIS) follows the expected shape, with two peaks during morning and evening rush hours reaching 10–12 and 9–11 ppbv in DJF and JJA, respectively. The diurnal cycle for the UCMF impact (DMET) is negative throughout the whole day, with the peak decrease during evening hours reaching −5 and −8 ppbv in DJF and JJA, respectively. In the case of the impact of modified dry deposition (DLU_D), it has a somewhat different pattern in the two seasons. In DJF, it is negative throughout the day with a strong peak during morning hours (−0.04 ppbv) and a smaller evening peak (−0.03 ppbv). In summer, this impact is positive almost during the whole day with two peaks during morning and early evening hours reaching 0.18–0.2 ppbv, while during night, the impact can be slightly negative up to −0.04 ppbv. The impact of BVOC changes (DBVOC) is very small during winter with negative values peaking at less than −0.01 ppbv. During summer, the impact is stronger with a clear positive peak during evening hours reaching 0.06 ppbv.

In the case of SO₂, the diurnal pattern for the impact of emissions and UCMF is similar to NO₂. The emissions impact peaks at morning and evening rush hours for JJA, reaching 2.6–2.8 ppbv, while in DJF, the maximum impact is reached at evening hours and the impact remains high during the whole night (around 2.5–3 ppbv). The DMET impact is negative with an evening peak reaching −0.06 and −0.03 in DJF and JJA, respectively. The impact on dry deposition is negative in JJA with maximum impacts during morning and evening hours reaching −0.05 to −0.07 ppbv. During JJA, the impact is positive during the day with increases of up to 0.03 ppbv and decreases during night of up to −0.04 ppbv. We already saw in the box plots and also expected that the impact of BVOC emission change has an almost zero effect on SO₂, which is not directly chemically tied to VOC chemistry.

Finally, for ozone, the impact of urban emissions is a decrease with two peaks during morning and evening hours reaching −10 to −12 ppbv in DJF and −8 to −10 ppbv in JJA. The impact of UCMF shows a clear increase peaking during evening hours, reaching around 5 and 10 ppbv during DJF and JJA, respectively. The impact of modifications of dry deposition is positive throughout the day with a strong peak during noon to early evening hours – during DJF, the peaks reach 0.2 ppbv, while a much stronger increase is modeled during summer reaching 1.5–2 ppbv. The impact of BVOC changes on ozone is virtually zero during DJF and negative during JJA, with a peak decrease around noon reaching 1 ppbv.

We also evaluated the diurnal cycle of the impact on deposition velocities, as this helps the interpretation of the DLU_D. In Fig. 11 the 2015–2016 winter and summer average of the DV diurnal cycle for the three analyzed pollutants is plotted as are the absolute values corresponding to the rural (“Nourban”) case. In the case of NO₂, DVs are reduced when turning rural land use into urban land use, and the maximum decrease occurs during noon to early afternoon, reaching −0.4 mm s⁻¹ in DJF with a stronger decrease reaching −3 mm s⁻¹ in JJA, while during night, the change is close to zero. Similar decreases are calculated for ozone with somewhat smaller nocturnal values in DJF and a bit weaker decrease during summer peak values. For SO₂, the impact on DV is different between DJF and JJA. During DJF, DV increases by 0.6 mm s⁻¹ during night, while a smaller increase is calculated around noontime (0.2 mm s⁻¹). During JJA, DV change for SO₂ is slightly above zero and a strong negative peak occurs during the day, reaching about −1.5 to −2 mm s⁻¹. Comparing with the absolute values, the impact of urban land use change in winter can reach −10 % to −20 % for ozone and NO₂, while it is 20 %–30 % for sulfur dioxide. In summer the decrease is even higher, reaching 50 % for ozone and NO₂, while for SO₂, the relative decrease is about 30 %–40 %.

Figure 11 The diurnal cycle of the rural deposition velocities (DV; solid lines; left y axis) and the impact of urbanization (dashed lines; right y axis) for 2015–2016 DJF (a) and JJA (b) for NO₂ (blue), SO₂ (brown), and O₃ (red) (mm s⁻¹).

4 Discussion and conclusions

We presented an analysis of the different contributors to the overall impact of urbanization (what we called here the rural–urban transformation; RUT) on urban gas-phase air pollutant concentrations. We focused on the four most important contributors to RUT, namely the impact of urban emissions (DEMIS), the impact of the urban canopy meteorological forcing on pollutant chemistry and transport (DMET), the impact of modified dry deposition due to the land cover modifications (DLU_D), and the impact of modified biogenic emissions due to modified land cover (and associated vegetation change) and modified meteorological conditions (DBVOC). By performing multiple simulations wherein each contributor of RUT was added one by one to the reference state representing land without urban land cover and urban emissions, we could quantify them individually.

The validation showed a reasonable range of model bias, and the annual cycles of pollutant concentrations for ozone and NO₂ were well captured. The same is true for the model ability to resolve the diurnal cycle of ozone. Regarding NO₂ biases, our results show a clear improvement from our previous study in Huszar et al. (2021), which used an almost identical setup and the same emissions input. It is clear that this improvement also cannot be explained by improved meteorology as the driving RegCM and CAMx simulations were the same. The only probable explanation is that the improvements were achieved by updating the chemistry mechanism in our simulations from CB5 to CB6r5. Indeed, CB6 was added to CAMx to take into account the long-lived organic compounds formed by peroxy radical reactions, which serve as an inhibitor of OH recycling and reduces NO_x removal by OH oxidation (Cao et al., 2021). Previously, Luecken et al. (2019) also found a better model performance for reactive nitrogen when CB6 was used instead of CB-V. The slight deviation of the monthly cycles from observed values is probably caused by incorrect annual temporal disaggregation profiles, which are dependent only on the emission activity sector but not on the geographic location. Some studies using older chemistry mechanism also obtained larger NO₂ biases (Karlický et al., 2017; Tucella et al., 2012), which suggests that the accuracy of the chemistry mechanism is probably very important. For ozone, monthly values were well represented by our model system, and the choice of the chemistry mechanism probably contributed to this, as in older studies using CB5 (or even the older CB-IV) (Zanis et al., 2011; Huszar et al., 2016a, 2020b) the biases were higher (moreover, the latter study used the same emission data as here) and were often caused by a strong nighttime bias, which seems to be partly removed in our study. Further, our results show similar model–observation agreement as the large online coupled model comparison study by Im et al. (2015). In the case of SO₂, the model is rather unable to correctly resolve the annual cycle of near-surface concentrations. We also saw this behavior in a similar manner in Huszar et al. (2016a) and in Karlický et al. (2017), and it points to deficiencies in the annual profile used to time-disaggregate annual emissions to monthly ones. The SO₂ biases can also be caused by wrong vertical turbulent mixing as large quantities of this pollutant are emitted from tall stacks and they have to mixed down to the surface layer, which is greatly influenced by the model representation of vertical eddy diffusivities. These are especially important in urban areas, and large uncertainty still persists in their calculation (Huszar et al., 2020a). In summary, for NO₂ and O₃ we did not identify substantial model biases in simulating urban near-surface concentrations of the analyzed pollutants. For sulfur dioxide, our model failed to correctly resolve the annual cycle, which suggests that the impact of urban SO₂ emissions can also be overestimated or underestimated depending on the model bias and should be perceived with caution.

The total impact of urbanization on NO₂ was calculated to around 3 (1 ÷ 5) ppbv in summer and 6 (3 ÷ 8) ppbv in winter. These numbers are smaller than the annual mean contributions calculated for 2001–2010 in Huszar et al. (2016a), and higher contributions were also modeled by Im and Kanakidou (2012); however, both simulated only the effect of urban emissions without considering the effect of the UCMF, which decreases near-surface concentrations (see further). The total impact on SO₂ is between 0 and 1 ppbv in summer and 0.5–1.3 in winter, which is a smaller contribution than in Huszar et al. (2016a) due to much lower sulfur emissions in 2015 compared to the 2005 emissions used there and due to not considering the UCMF effects in the earlier study. The total average contribution for ozone summer MDA8 is about 1.5 (0 ÷ 2) ppbv. In Huszar et al. (2016a) for an ensemble of central European cities and Im et al. (2011a, b) for Mediterranean cities a decrease in ozone was shown (and increase over rural areas, similar to our results), but they accounted for only the urban emission impact. Indeed, urbanization via the UCMF increases ozone concentrations (Kim et al., 2015; Huszar et al., 2018a, 2020b), which can offset the decrease seen solely due to urban emissions. For all three pollutants, the effect of emissions (DEMIS) is stronger than the total effect of urbanization (DTOT) due to the strong modulating effect of the urban canopy meteorological forcing. As already calculated by many (e.g., Wang et al., 2007, 2009; Struzewska and Kaminski, 2012; Zhu et al., 2015; Huszar et al., 2020a), vertical eddy diffusion is the most important component of UCMF, which is strongly enhanced above urban areas. Consequently, it leads to reduced near-surface concentrations of primary pollutants (e.g., NO₂, SO₂) and an increase in ozone due to reducing the titration by NO (Escudero et al., 2014; Xie et al., 2016a, b).

As for the impact of UCMF (DMET) alone, our simulations showed a decrease by about 2 ppbv for NO₂ and by 0.2–0.3 ppbv for SO₂ (for both seasons) as well as an increase in summer MDA8 ozone by about 2 ppbv. These numbers fit previous findings in Sarrat et al. (2006), Struzewska and Kaminski (2012), Kim et al. (2015), and Huszar et al. (2018a, 2020b) well. They concluded that three main components play the most important role in UCMF: increased urban temperatures, decreased wind speeds, and increased vertical turbulent diffusion. While elevated surface temperatures favor photochemistry, they also result in stronger dry deposition as shown by Huszar et al. (2018a). Regarding the wind speed and turbulence effect, they are counteracting, which is seen in our results too, especially for SO₂. For some of the cities, the impact is positive, meaning that the reduction of winds results in the emitted material remaining close to the sources. This was also previously seen by Huszar et al. (2018b) wherein the turbulence and wind effects were strongly competing. Our results also showed that the trade-off between wind and turbulence effects also depends on how the model simulates the UCMF components, and in our results, WRF produced a somewhat stronger increase in turbulence due to UCMF and weak wind reduction compared to RegCM. For ozone, the UCMF increased ozone by 2 ppbv, which is in line with previous finding in Huszar et al. (2018a), although they also included the effect of BVOC emissions modifications, which was treated separately here (see further). Due to urbanization, a similar increase was obtained by Martilli et al. (2003), Jiang et al. (2008), Xie et al. (2016a), and Jacobson et al. (2015). Some authors found somewhat larger increases for ozone (e.g., Ryu et al., 2013, for Seoul), but they adopted higher resolutions for the cities in focus and thus obtained higher peak impacts in urban centers (as seen in, e.g., Huszar et al., 2020b, too).

The diurnal cycle for the DMET impact shows a very characteristic pattern. In the case of primary pollutants (NO₂ and SO₂) the decrease is strongest during evening hours. This can be explained by the largest absolute values during evening hours, which is further closely related to the maximum impact of emissions – a similar finding was found by Huszar et al. (2018a) and also by Huszar et al. (2018b) for primary aerosol components. Indeed, the quantity of gases transported due to enhanced turbulence is proportional to the absolute concentrations, and these are highest during evening hours due to strong emissions during transport rush hours. For ozone, the diurnal pattern contains a maximum during evening hours corresponding to the largest impact on NO₂. This justifies the argument that the UCMF-induced ozone increase is mainly caused by reduced NO_x due to strong urban dilution and consequent reduced titration.

Besides the strong and well-documented air quality effects of urban emissions (DEMIS) and UCMF (DMET), our study also looked at two other contributors to RUT, which were expected to be smaller but which were not yet quantified in detail. Our study, at least according to the knowledge of the authors, is among the first to explicitly investigate the effect of urbanization from the perspective of change in dry deposition (DLU_D), and we also looked at the effect of urbanization-induced changes in BVOC emissions, which was examined only partly in previous studies (e.g., Huszar et al., 2018a; Li et al., 2019).

The impact due to modified dry deposition shows a distinct picture for NO₂ between summer and winter. For both seasons, reduced deposition velocities were modeled with stronger decreases in summer (the WRF-driven CAMx results are not shown as they differ from the RegCM-driven only slightly). Reduced deposition velocities result in higher concentrations, which is opposite to what was modeled. To better understand what controls the NO₂ budget, we have to consider the simultaneous effect of ozone changes due to changes in dry deposition. Our results showed strong increases in ozone concentrations, which is probably caused by suppressed dry deposition (for winter too; not shown in this paper). This is expected as many studies showed strong dependence of ozone concentrations on ozone deposition (Tao et al., 2013; Park et al., 2014) and ozone deposition on land use information (Mcdonald-Buller et al., 2001). When examining the concentration response to changed dry deposition, one has to consider the indirect impact due to other influenced pollutants, and pollutants responsible for NO₂ removal (e.g., by reaction NO₂+OH, forming nitric acid, or by NO₂+O₃, forming the nitrate radical) were probably impacted by weaker dry deposition (as seen for ozone), resulting in decreases in NO₂, outweighing the direct impact of dry deposition change, as seen for winter. Another factor playing a role in decreases in NO₂ can be the much larger (by 50 %) deposition velocities for nitric acid (HNO₃) in the Zhang model for the urban land use type compared to crops or similar rural land use (i.e., “Nourban” case). Large dry deposition for HNO₃ in turn results in a decrease in this compound, which reduces the recycling of NO₂ from it (by photolysis). On the other hand, in summer, such effects can amplify the impact. In this season the deposition-induced ozone changes probably played a role in the NO₂ budget. It has to realized that a major pathway of NO₂ chemistry in cities is oxidation of NO with ozone ($\mathrm{NO}+{\mathrm{O}}_{\mathrm{3}}\to {\mathrm{NO}}_{\mathrm{2}}$). Increased ozone concentrations thus result in more NO oxidizing to NO₂.

In the case of sulfur dioxide, deposition increased in winter, which resulted in a clear decrease in near-surface concentration. The dry deposition of sulfur strongly differs between wet and dry soils (Hardacre et al., 2021), and according to Zhang et al. (2003), who provided the dry deposition scheme we adopted, the deposition velocities (DVs) are higher for urban areas than for crops or similar rural land use types (which was considered in the “Nourban” case). In winter, soils are very often wet, which could result in an increase in DVs (especially during night as seen in our results) and a consequent decrease in concentrations. During summer, DVs decreased for SO₂; however, there is no clear increase in concentrations, i.e., almost no change in RegCM-driven simulations and even some decrease in the WRF-driven ones. This again can be explained by the impact of deposition on other chemical species which cause removal of SO₂, typically oxidation by the OH radical.

The impact of BVOC emission changes (DBVOC) is straightforward and expected for ozone, i.e., a decrease by 0.5–1 ppbv. BVOC emissions decreased due to urbanization-related reduction of vegetation (i.e., reduction of vegetation fraction and leaf area index) and increased due to higher urban temperatures (within the action of the UCMF). This latter effect was smaller, resulting in the dominance of the first effect and an overall decrease in emissions, which is a similar result as in, e.g., Li et al. (2019). As ozone chemistry in cities in Europe (and also North American and Asian megacities) is characterized by a VOC-controlled regime with a high ${\mathrm{NO}}_x/\mathrm{VOC}$ ratio (Beekmann and Vautard, 2010; Xue et al., 2014), ozone quickly responds to changes in VOC emissions; i.e., it decreases with decreasing BVOC emissions. This is in accordance with previous studies: e.g., Song et al. (2008) reported an 10 % decrease in ozone concentrations. The reduction in ozone was shown to be largest during daytime, which is in accordance with the largest BVOC emissions. Previously, Huszar et al. (2018a) reported ozone increases due to BVOC changes due to UCMF alone (i.e., not considering the impact of reduced vegetation) of the order of up to 0.1 ppbv. Our study showed that if vegetation modifications related to urbanization are added, this increase is outweighed by a much stronger decrease due to lower BVOC emissions.

Simultaneously with the DBVOC-related ozone decrease, we calculated a small summer increase in NO₂ by about 0.01 ppbv. This cannot be explained by potentially reduced NO concentrations and suppressed NO₂ formation with the reaction of ozone (titration), as NO also increased slightly as the result of DBVOC (not shown explicitly in this study). Moreover, soil NO_x emissions in MEGAN also decreased slightly due to urban land use transformation. Reduced BVOC emissions result in reduced peroxy radical (RO₂) concentrations, which is an important oxidation pathway to form NO₂ from NO (Geng et al., 2011) and would result in a decrease in NO₂. There must therefore be another compensating mechanism responsible for NO_x increase, and this is probably the reduced concentrations of NO_x sinks. One of the important urban contributors to this are PANs (peroxy acetyl nitrates), and as biogenic VOCs are a major contributor to urban PAN concentrations, it can be expected that with decreased BVOC emissions, the PAN sink is reduced, resulting in higher NO_x concentrations (Fischer et al., 2014; Toma et al., 2019). Another reason can lie in reaction with the OH radical, which is reduced if ozone is reduced. In short, the relatively small positive NO₂ response to urbanization-induced biogenic emissions changes is probably a simultaneous action of multiple indirect chemical pathways, and deeper process-based analysis should be performed to explicitly show the contribution and trade-off of each of them. The strongest changes modeled for Milan, Italy, can be explained by its relatively warm climate and large size, making the BVOC emission reduction strong and thus having a stronger effect on ozone and NO₂ (via secondary effects discussed above).

The diurnal patterns of different RUT contributors are explainable by the typical diurnal patterns of anthropogenic emissions and also by the UCMF's diurnal cycle. For the impact of urban emissions only, the DEMIS impact clearly resembles the emissions which peak in cities at morning and evening rush hours (Huszar et al., 2021). The impact of urbanization-induced meteorological changes (UCMF) is governed by increased vertical eddy diffusion, which is strongest during daytime (Huszar et al., 2020a) and causes a decrease in primary pollutant concentrations as well as an increase in ozone (reduced titration; Huszar et al., 2018a). As the impact is proportional to the absolute concentrations, it is again largest during morning and (especially) early evening rush hours, during which even stronger vertical mixing usually occurs. A somewhat more complex explanation is needed regarding the diurnal pattern of the land use change and dry deposition impact. For ozone, a strong decrease occurs during daytime, which is related to high daytime absolute values and a stronger decrease in deposition velocities during the day due to urbanized land surface. Over vegetated surface, ozone dry deposition is determined mainly by stomatal resistance, which is lowest during daytime, meaning that the ozone dry deposition velocities are highest during the day (Park et al., 2014). When the surface is urbanized, the vegetation's role in dry deposition becomes very small and this daytime peak of DV disappears, which results in the strongest decrease in DV during the day. The diurnal cycle of the DLU_D for NO₂ is probably the combined effect of decreased DVs, which are strongest during the day and logically lead to increased concentrations, and the reduction due to a dry-deposition-induced increase in ozone concentrations (see above regarding how this might influence the NO₂ concentrations). As this latter occurs during noontime, the resulting shape of the diurnal pattern of DLU_D for nitrogen dioxide has a double peak. Finally, the diurnal cycle of the response of SO₂ to dry deposition changes due to urban land surface has a clear daytime peak during winter when the impact is smallest. This is in line with the smallest increase in DV for SO₂, which occurs during daytime. During summer, the shape of the diurnal cycle for the impact on SO₂ is similar; however, the values are higher, making the peaks reach positive values. These can be attributed to the summer decrease in DV, peaking during daytime, which leads to an increase in concentrations. A more interesting feature is that the impact goes to negative values during nighttime for this pollutant. This cannot be explained by the DVs alone, and, as already said above, some secondary chemical effects play a role; e.g., increased ozone concentration leads to increased OH radical, which can lead to increased oxidation of SO₂, as already suggested earlier. The diurnal impact of DBVOC for ozone in summer has a clear daytime negative peak, which is connected to the highest BVOC emissions occurring at this time of the day. As the most important factor in the urbanization impact on BVOC is the reduction of vegetation, the highest reduction of emissions is expected to occur during daytime, and in a VOC-limited environment this results in a stronger suppression of ozone production (i.e., leading to reduction). In the case of NO₂, the DBVOC impact is probably connected to decreased PAN, which is an important sink for NO_x (as already said above). As the absolute NO_x concentrations peak during evening rush hours, this sink has the strongest effect during this time, creating the evening peak of DBVOC for NO₂.

Our results further showed that the magnitude of each contributor for Prague is higher at the 1 km resolution nested domain than at 9 km resolution. This was already seen in previous studies focusing on the impact of UCMF in extreme air pollution (Huszar et al., 2020b) and also concerned the magnitude of the UCMF itself (e.g., the temperature increases due to urbanization have a higher city core peak at higher resolution). This was further seen regarding the impact of urban emissions only in Huszar et al. (2021). It can be assumed that if a nested domain approach at high resolution was applied to other cities, the impacts would have a stronger city core peak than at 9 km resolution.

To summarize our finding, we showed with an ensemble of 19 central European cities that the strongest contributors to the impact of rural-to-urban transformation are the urban emissions themselves (increase concentrations for nitrogen dioxide and sulfur dioxide and decrease those for ozone) and the urban canopy meteorological forcing (decreases the concentration of primary pollutants and increases those of ozone). They have to be accounted for simultaneously as the impact of urban emissions without considering UCMF can lead to overestimation of the impact (Huszar et al., 2021). Additionally, we quantified two weaker contributors. The effect of modified land use on dry deposition and the effect of modified biogenic emissions have 1 order weaker magnitude than emissions and the UCMF. However, we showed that for summer ozone, these are strong and of comparable order as the two strongest impacts. In other words, when analyzing the overall impact of urbanization on air pollution for ozone, all four contributors have to be accounted for, having a similar order of magnitude, while for primary gas-phase pollutants (i.e., NO₂ and SO₂), the two weaker contributors are by at least 1 order of magnitude smaller and the error made is small if they are neglected.

Finally, we must stress that we focused on cities from a relatively small region, meaning the cities do not constitute substantially different climates and the typical “rural” vegetation was considered to be crops. In other parts of the world, e.g., tropical areas, rural-to-urban transformation takes place over different vegetation cover and, e.g., the impact of modified BVOC emissions could be much stronger than in the case of central European cities. Further, it has to be noted that some secondary effects of modified pollutant concentrations can potentially also play a role via the direct and indirect radiative effect of emissions. The direct effect of aerosols can alter photolysis rates and temperatures, influencing air chemistry (Han et al., 2020; Wang et al., 2022). It was further shown by many that aerosol emitted by urban areas modulates the vertical structure of the atmosphere, resulting in modification of stability and/or convection (Miao et al., 2020; Slater et al., 2022; Fan et al., 2020; Yu et al., 2020), which in turn can modify the vertical mixing or precipitation (Zhou et al., 2020; López-Romero et al., 2021); this finally feeds back to influence species concentration via wet deposition and mixing. Our study was an offline coupled one, meaning that no feedbacks from species concentrations via radiation and cloud–rain microphysics were accounted for. These studies, however, indicate that to obtain an even more comprehensive picture of the total RUT impact, these secondary effects also have to be taken into account in the future.