Development of an atmospheric chemistry model coupled to the PALM model system 6.0: Implementation and ﬁrst applications

. In this article we describe the implementation of an online-coupled gas-phase chemistry model in the turbulence resolving PALM model system 6.0. The new chemistry model is part of the PALM-4U components (read: PALM for you; PALM for urban applications) which are designed for application of PALM model in the urban environment (Maronga et al., 2020). The latest version of the Kinetic PreProcessor (KPP, 2.2.3), has been utilised for the numerical integration of gas-phase chemical 5 reactions. A number of tropospheric gas-phase chemistry mechanisms of different complexity have been implemented ranging from the photostationary state to more complex mechanisms such as CBM4, which includes major pollutants namely O 3 , NO, NO 2 , CO, a simpliﬁed VOC chemistry and a small number of products. Further mechanisms can also be easily added by the user. In this work, we provide a detailed description of the chemistry model, its structure along with its various features, input requirements, its application and limitations. A case study is presented to demonstrate the application of the new chemistry 10 model in the urban environment. The computation domain of the case study is comprised of part of Berlin, Germany, covering an area of 6.71 x 6.71 km with a horizontal resolution of 10 m. We used "PARAMETERIZED" emission mode of the chemistry model that only considers emissions from trafﬁc sources. Three chemical mechanisms of varying complexity and one no-reaction (passive) case have been applied and results are compared with observations from two permanent air quality stations in Berlin that fall within the computation domain. The results show importance of online photochemistry and dispersion of air 15 pollutants in the urban boundary layer. The simulated NO x and O 3 species show reasonable agreement with observations. The agreement is better during midday and poorest during the evening transition hours and at night. CBM4 and SMOG mechanisms show better agreement with observations than the steady state PHSTAT mechanism.

. Schematic of PALM-4U chemistry model Depending on the need, a user can select a chemistry mechanisms of different complexity. The Fortran code for the selected gas-phase chemistry mechanism is generated by a preprocessor based on KPP (Damian et al., 2002). The latter is described in more detail in Section 2.2.2. Besides chemical transformations in the gas phase and a simple photolysis parameterization (Section 2.2.3) the chemistry module includes dry deposition (Section 2.2.5), an interface to the aerosol module SALSA (Kurppa et al., 2019) (Section 2.2.4) and an option for anthropogenic emissions (Section 2.3). 5 5 https://doi.org /10.5194/gmd-2020-286 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.

Prognostic equations
When gas phase chemistry is invoked, N additional prognostic equations are solved, with N being the number of variable compounds of the chemical reaction scheme. Except for the SGS flux terms, the overbar indicating filtered quantities is omitted to improve readability. The three-dimensional prognostic equation for an atmospheric pollutant then read as: 5 where c n (n = 1, N ) is the concentration of the respective air constituent, which can be either a reactive or passive gas-phase species or an aerosol particulate matter compound. The term on the left-hand side is the total time derivative of the pollutant concentration. The first two terms on the right-hand side represent the explicitly resolved and the SGS transport of the scalar chemical quantity in x, y and z direction. A double prime indicates a SGS variable. The third term represents the change in concentration (c n ) of the trace gas n over time due to production and loss to chemical reactions, which can be described as (VOCs) and a very small number of products. For the convenience of the users, it is not required to run kpp4palm for the ready-to-use mechanisms (Table 1) as their Fortran subroutines are already supplied with PALM. Currently PHSTATP is the default mechanism which will automatically be compiled with the rest of the PALM source code when the chemistry option is switched on. However, users can also add modified versions of the existing chemical mechanisms or define completely new mechanisms according to their specific needs.

Photolysis frequencies
The parameterization of the photolysis frequencies is adopted from the Master Chemical Mechanism (MCM) v3 according to Saunders et al. (2003). Photolysis frequencies are described as a function of the solar zenith angle ϑ and three parameters, which are specific for each photolysis reaction: Values for l, m and n are given for the relevant photolysis reactions in Saunders et al. (2003) and on the MCM web page (http://mcm.leeds.ac.uk/MCM/parameters/photolysis_param.htt). Currently only a simple parameterized photolysis scheme is available for photochemical reactions. More extensive photolysis schemes such as Fast-J based on radiative transfer modelling will be included in the future. These models will also make use of the shading due to buildings, which is already implemented for the shortwave radiation in the PALM-4U urban surface model.

Coupling to SALSA aerosol module
The sectional aerosol module SALSA2.0 (Kokkola et al., 2008) has recently been implemented into PALM and a detailed description is given in Kurppa et al. (2019). SALSA describes the aerosol number size distribution, aerosol chemical compo-

Deposition
Deposition is a major sink of atmospheric pollutant concentrations. Currently, only dry deposition processes are included as precipitation (leading to wet deposition of pollutants) is not yet included in PALM. For dry deposition, a resistance approach is taken where the exchange flux is the result of a concentration difference between atmosphere and earth surface and the resistance between them. Several pathways exist for this flux, each with its own resistance and concentration. The aerodynamic 5 resistance depends mainly on the atmospheric stability. In PALM, it is calculated via Monin-Obukhov similarity theory, based on roughness lengths for heat and momentum and the assumption of a constant flux layer between the surface and the first grid level.
For gases, the quasi-laminar layer resistance depends on the atmospheric conditions and diffusivity of the deposited gas and it is calculated following (Simpson et al., 2003). Finally, the surface (canopy) resistance for gases, which is the most 10 challenging resistance to estimate due to the enormous diversity of surfaces, is calculated using the DEPAC module (Van Zanten et al., 2010). DEPAC is widely used in the flux modelling community (e.g. (Manders et al., 2017;Sauter et al., 2018).
The surface resistance parameterizations are different for different land use types defined in the model. In DEPAC, three deposition pathways for the surface resistance are taken into account: through the stomata 15 through the external leaf surface through the soil DEPAC is extensively described in a technical report by Van Zanten et al. (2010). It also includes a compensation point for ammonia which is currently set to zero in PALM.
For the passive particulate matter in the chemistry model, the land-use dependent deposition scheme of Zhang et al. (2001) 20 has been implemented into PALM. The formulations have been chosen as they include an explicit dependence on aerosol size. For particulate matter, the deposition velocity is calculated by the gravitational settling or sedimentation velocity (mainly relevant for the larger particles), the aerodynamic resistance (see above) and the surface resistance. The sedimentation velocity mainly depends on particle properties, the gravitational acceleration and the viscosity coefficient of air. The formulation for the surface resistance is empirical with parameters that are based on a few field studies including the collection efficiencies for 25 Brownian diffusion, impaction and interception, respectively, and a correction factor representing the fraction of particles that stick to the surface depending on the surface wetness. Further details can be found in Zhang et al. (2001).

Traffic emissions
The chemistry model of PALM includes a module for reading gaseous and passive anthropogenic emission input from traffic sources and converting it to the appropriate format. These emission data can be provided in three possible levels of detail 30 (LODs), depending on the amount of information available at the user's disposal. With LOD 0, traffic emissions are parameterized based on emission factors specific to particular street types. All street segments contained in the domain will be classified into "main" and "side" street segments. A mean surface emission flux tendency for each chemical species contained the active chemical mechanism, in kilograms per square meter per day, will be provided together with a weighting factor for main and side street emissions in the PALM parameter file. Street type classification based on Open StreetMap definitions (Open-StreetMap contributors, 2017) is to be included in the PALM static driver (Maronga et al., 2020). A diurnal profile derived from traffic counts is implemented to disaggregate total emissions into hourly intervals. Currently a default profile is applied to 5 all species for main and side street segments, and details can be found in the on-line documentation for the PALM-4U chemistry model (https://palm.muk.uni-hannover.de/trac/wiki/doc). More detailed traffic emission data can be provided in gridded form in PALM-specific netCDF format (Maronga et al., 2020). LOD 1 accepts a spatially distributed annual emission data for each discretized street cell in the computation domain, which will be temporally disaggregated using sector-specific standard time factors. In addition, using LOD 2, the user can further introduce preprocessed emission data that are already temporally 10 disaggregated, in hourly intervals, for instance. Future plans include expansion of the emission model to accommodate further modes of anthropogenic emissions such as domestic heating.

Initial and boundary conditions
Lateral boundary conditions for chemical compounds can be chosen in the same way as the lateral boundary condition for other scalars, e.g. potential temperature, being either cyclic conditions or non-cyclic (Maronga et al., 2020). In most urban 15 applications, chemistry requires non-cyclic boundary conditions, because cyclic conditions lead to accumulation of pollutants to the modelling domain if pollutant emissions exist. As part of the PALM-4U components, nesting has been implemented to chemistry module. In offline nesting, PALM can be coupled to a larger (meso, regional, or global) scale model to provide dynamic boundary conditions for the meteorological variables as well as air pollutants. As larger scale models do not fully resolve turbulence, a synthetic turbulence inflow generator has been introduced (Gronemeier et al., 2015).

20
Initial concentrations of primary compounds to control the chemistry model are controlled by the chemistry namelist "&chemistry_parameters" in PALM parameter file. Options for prescribing initial conditions are available for both surface initial conditions and initial vertical profiles for the area or region of interest. There are no default initial concentrations and the user is responsible for providing these values based on, e.g., measured background concentration of primary chemical compounds for whom initial concentration are defined. These primary compounds must be part of the applied chemical mechanism. 25 3 Chemistry model application In order to demonstrate the ability of the chemistry implementation, we performed a simulation for an entire daily cycle in a realistic urban environment and compared simulation results against observational data. To analyse the effect of different chemical mechanisms with different complexity on the resulting concentrations, we performed simulations with three chemical mechanisms and one passive case where only transport and dry-deposition was considered. 30 10 https://doi.org/10.5194/gmd-2020-286 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.

Numerical setup description
The PALM model system, version 6.0, revision 4450 and 4601 (only for flux profiles of chemical compounds), has been used in this study. Details of the dynamic core of the model are described in Maronga et al. (2015) and Maronga et al. (2020). The model was run for 24 hours from 00:00 UTC (0200 CEST 17 July 2017), to 24:00 UTC (02:00 CEST 18 July 2017) for a city quarter located around the Ernst-Reuter-Platz in Berlin, Germany. This particular day was chosen as it represents an 'ideal' 5 Berlin summer day with mostly sunny with scattered clouds in the morning and passing clouds in the afternoon and evening hours. The temperature ranged between 289 K to 298 K with moderate winds predominantly from westerly direction. It was not a weekend, therefore, emissions from the traffic were not affected by the reduced traffic. In addition to near surface observations from online web sources Thorsen (2017), we also analysed radiosonde data from Lindenberg (

15
Observations from radiosonde and Ceilometere showed strong showers on the previous days that resulted in a very well mixed (almost moist adiabatic) layer in the lower troposphere that lead to an almost constant potential temperature gradient above the inversion, therefore, there was no residual layer at mid-night on 17 July, when the model was initialised for a 24 hour run. The Ceilometer observations for 17 July, 2017, also did not show disturbances by clouds, the mixed layer top however, remained below 2000 m throughout the diurnal cycle. Fig. 2 shows the computational domain that covers an area of 6.71 km 20 x 6.71 km (671 x 671 grid points) with a model top at 3.6 km above surface. The horizontal grid spacing in x and y direction is 10 m. In the vertical, with 312 layers, the grid spacing is 10 meters up to 2700 m, above which it increases by an expansion factor of 1.033 until the grid spacing in the z-direction reaches 40.0 m.
Rayleigh damping has been used above 2500 m in order to weaken the effect of gravity waves above the boundary-layer top.
A multigrid scheme has been used to calculate pressure perturbations in the prognostic equations for momentum (Maronga 25 et al., 2015). A third order Runge-Kutta scheme (Williamson, 1980) has been used for time integration. The advection of momentum and scalars was discretized by a fifth-order advection scheme by Wicker and Skamarock (Wicker and Skamarock, 2002). Following Skamarock (2006); Skamarock and Klemp (2008), we employed a monotonic limiter for the advection of chemical species along the vertical direction in order to avoid unrealistically high concentrations within the poorly resolved cavities (e.g. courtyards represented by only a few grid points) which can occasionally occur due to stationary numerical   The profiles of potential temperature indicate a vertically well-mixed boundary layer during daytime, evolving from approximately 500 meters at 7:00 CEST to more than 2300 m at 21:00 CEST, while in the evening hours the near-surface layer stabilises. The mixed-layer depth agrees fairly well with the observed values from Ceilometer measurements (horizontal bars 30 in Fig. 3a), except for the late afternoon and evening hours, where the modelled boundary-layer depth is over-predicted by Throughout the diurnal course, near surface concentrations of NO are 1 to 2 ppb higher than the remaining mixed layer. At 9:00 CEST owing to emissions from traffic sources, the NO concentration increased from 3 ppb to 5 ppb in the first few meters above surface, then evenly dilutes in the shallow mixed layer and stabilises with a mixing ratio of 3 ppb. Besides the high NO x emissions in the morning hours this increase can also be attributed to the onset of NO 2 photolysis (R2). During rest of the day,  Unlike NO, the NO 2 concentration profiles in Fig. 4b show only small differences in the morning hours (7:00 and 9:00 CEST). In the afternoon (12:00 and 17:00 CEST) when the convection is stronger, the NO 2 concentration is the lowest due to the combined effect of upward vertical mixing indicated by Fig. 4f and photolysis of NO 2 (R2).
In addition to reaction R3, the O 3 levels within the mixed layer are also replenished through down-welling from above the 10 inversion during the day which is evident from the negative flux profiles of O 3 Fig. 4g. As a result O 3 concentration gradually increased from 18 ppb at 07:00 CEST to 43 ppb at 17:00.
In the evening hours (21:00 hour on wards), NO concentration reduced to 1 ppb near the surface while it is completely removed from the residual layer. In the reduced(no) solar radiation, NO 2 photolysis (R2) slows-down(stops) whereas emission and limited R1 reaction (because of very low available NO) favours increase in the NO 2 mixing ratio. This resulted in the highest near surface NO 2 concentrations at 24:00 CEST in the nocturnal stable layer. In the well mixed residual layer above, the NO 2 concentration at 24:00 CEST is up to 1 ppb lower than NO 2 concentrations at 21:00 CEST. Ozone also reduced up to 3 ppb at 24:00 CEST but remained well mixed in the residual layer above the shallow nocturnal stable boundary layer. Carbon monoxide has been added to the analysis as a proxy to passive species Fig. 4d . In the morning hours concentration is the 5 highest, that gradually reduced due to turbulent mixing in the rapidly growing mixed layer and due partly to dry deposition.

Spatial distribution of pollutants
The spatial distribution of NO 2 and O 3 concentrations 5 m above surface (Fig. 6) are discussed for 9:00 and 21:00 CEST as simulated with the CBM4 mechanism. The concentrations of NO x , aldehyde and hydrocarbon species (not shown), which are emitted by road traffic, were high in the morning and evening hours due to stable conditions and high emissions.

15
Most of the high NO 2 concentrations in the morning hours (Fig. 6a) are predicted in the street canyons in the southern part of the simulation domain, where numerous buildings with a height of approximately 20 m are present (see Fig. 2a). In the street the street canyons is inhibited mainly due to delayed heating of the ground which is attributed to the shading of the surrounding buildings. Therefore O 3 is predominantly titrated by reaction R1 (Fig. 6b)  In the evening, however, NO 2 distribution is somewhat more uniform over the entire road network (Fig. 6c). This is attributed to emissions from traffic with reduced or no photolysis of NO 2 (R2) and titration of O 3 . Under the low wind conditions below 50 m, ventilation is reduced and leads to increase in NO 2 over street crossings, main and side roads. Consequently, the NO 2 concentration increased by 30% in the evening hours over roads and adjacent paved areas whereas over the vegetated areas (grass, crops, shrubs and trees), NO 2 concentrations is around 5 ppb in To provide an overview of the pollutant dispersion at the street level, a small section of the model domain (500 x 500 m) in the Hardenbergplatz area has been analysed (Fig. 7). The brown filled dot shows the location of Hardenbergplatz air quality station. This small urban section is characterised by typical urban environment with streets, paved areas, buildings with varying heights, and some vegetated areas of the Tiergarten located to the north-east. The dispersion of air pollutants in a street 20 canyon generally depends on the aspect ratio (building height to street width ratio) and rate at which the street exchanges air vertically with the above roof-level atmosphere and laterally with connecting streets (N'Riain et al., 1998). Figure    Due to the photostationary state, there is no net O 3 production in PHSTAT mechanism ( Fig. 9a and b). In the morning hours despite reaction R1 that consumed most of the NO, the NO 2 concentration is slightly reduced to 7 ppb by 7:00 CEST compared to its initial near surface background levels of 10 ppb. The small drop in NO 2 concentration is attributed to the mixing of NO 2 within the shallow mixed layer reaching 500 m in depth. Titration of O 3 by NO, consequently reduces O 3 concentrations to 15 ppb. After 07:00 CEST, photodissociation leads to decrease of NO 2 to very low levels (2-4 ppb) 50 m above surface. However, 10 NO 2 concentration remains on the order of 6 to 7 ppb close to the surface due to emissions from roads and NO x -VOC-HO x chemistry. Consequently O 3 concentration peaks between 13:00 to 18:00 CEST by reaction R3. Above the rapidly growing mixed layer, the background concentration remained 2 and 40 ppb for NO 2 and O 3 respectively. In the evening and night time emissions from traffic, suppression of NO 2 photolysis (R2) and a slow R1 reaction increases NO 2 concentrations near the surface up to 20 ppb, while O 3 concentration is suppressed to less than 15 ppb due primarily to reaction R1. Above the stable 15 nocturnal boundary layer the NO 2 is well mixed in the near neutral residual layer with a concentration of around 10 ppb, while due to absence of NO in the residual layer, and downward transport of background O3 concentration through the entrainment zone, the O 3 concentration increases by 30 ppb in the residual layer.
Photolysis of NO 2 , daytime NO x -VOC chemistry in the case of the SMOG (Fig. 9c and d) and CBM4 ( Fig. 9e and  In the CBM4 mechanism NO 2 levels are reduced as much as 40% between 12:00 to 19:00 CEST whereas O 3 concentration increases with almost the same percentage. Since no reactions leading to a net production of O 3 are considered for the photostationary state (PHSTAT mechanism), the maximum ozone concentrations is determined by the photolysis of NO 2 and is 25 therefore lower than for the SMOG and CBM4 mechanism.
In all three mechanisms, a gradual decay of O 3 is found in the evening (after 20:00 CEST). Near the surface, the interaction of NO x with O 3 chemistry is more pronounced. During the daytime specially in the morning hours , down-welling of O 3 in the entrainment zone also contributes to increase in O 3 concentration in the mixed layer. In the evening, after the collapse of the mixed layer, a well-mixed residual layer with near-uniform pollutant concentration can be found above the shallow 30 stably stratified nocturnal boundary layer. Although behaviour of NO x and O 3 is essentially the same in all three mechanisms, however, the difference is more evident in CBM4 because we see a more pronounced diurnal course of O 3 in CBM4.

Comparison of pollutant concentrations with observations
The observation data from Wedding and Hardenbergplatz is compared with the mean concentrations of the modelled NO, NO 2 , and O 3 from five points in the vicinity of the Wedding and Hardenbergplatz stations (Fig. 10). Error bars indicate standard deviation at every hour in the time series data from the five points. Ozone is not measured at Hardenbergplatz, and therefore, the time series curve of the observed O 3 is not included in Figure 10 and the analysis at Hardenbergplatz is restricted to NO 5 and NO 2 species only.
The three mechanisms (PHSTAT, SMOG and CBM4) are able to reproduce the diurnal cycle of NO, NO 2 , and O 3 with reasonable accuracy. CBM4 shows a more pronounced diurnal course than SMOG and PHSTAT. The timing of the maximum ozone concentration is reproduced much better in CBM4 than SMOG and PHSTAT. Since PHSTAT does not include any net O 3 formation, the maximum O 3 concentration occurs at the time of maximum solar elevation. However, SMOG and CBM4 10 include photochemical ozone production which leads to higher maximum ozone concentrations. A shift of the maximum ozone concentration from local noon to later afternoon is in agreement to observations. Due to a more detailed description of NO x -VOC-HO x chemistry this process is described more realistically for CBM4 than for SMOG. All three mechanisms failed to realistically simulate evening peaks in NO x (NO and NO 2 ) concentrations and the corresponding low O 3 concentrations. This can possibly be related to the slow cooling of the surface as described earlier which results in less pronounced stable layer near 15 the ground after sunset. Since O 3 is inversely related to NO x , the underestimated NO x -concentration in the evening leads to reduced titration of O 3 , therefore, results in higher levels of O 3 . Another important reason for bias in the model is the missing emission sources. We utilized emission information only from traffic sources, however, contribution from industry, household and biogenic volatile organic compounds (BVOC) is ignored which is a main source of uncertainty in the model.
The no-reaction case only contains transport and dry deposition of NO x . Without chemical reactions, NO continues to 20 increase, which is primarily due to the lack of chemical sinks in combination with the application of cyclic boundary conditions.
In contrast to the cases with chemical reactions, NO 2 is lower than NO since there is no conversion of NO to NO 2 . Despite titration by NO, a very slow increase in the near surface O 3 in the morning hours can be attributed to the downward mixing from the residual layer during the growth of the mixed layer. because reaction rates of chemical species vary greatly that leads to a time-step much smaller than meteorological parameters, therefore the computation speed of the entire system is largely effected. The simulation of the chemical reactions can take as much as 90% of the total computational time in the calculation of an online-coupled simulation (Cao et al., 2019). The computation time of chemistry increases with increasing number of the chemical compounds and the number of chemical reactions in a given chemical mechanism. Table 2 shows comparison of the four chemical mechanisms relative to the meteorology only resource requirement of chemistry coupled LES simulations, it is important to put extra effort in the design of the simulation and selection of the appropriate chemical mechanisms.

Concluding remarks
We have outlined the structure and important features of a chemistry model that has been added to the PALM model system 6.0 as part of PALM-4U (PALM for urban applications) components and coupled to the PALM model core. PALM model system  The results of the presented case study demonstrate the ability of chemistry model to predict photochemical reactions, along with advection and deposition of pollutants. The difference between reactive and non-reactive case clearly indicates that invoking reactive chemistry is of critical importance to accurately predict pollutant concentrations at the micro-scale.
The simulated concentrations follow the observed diurnal cycle of the pollutants. However, the weak agreement warrants an improvement in the description of the parameterized emissions, the large-scale forcing and model resolution.

15
Especially with coarse grid resolution some street canyons are only poorly resolved by the numerical grid, occasionally leading to the situation that street canyons are only one grid point wide or isolated cavities occur. However, within such enclosed cavities the flow is only poorly resolved too, meaning that the total transport of energy and matter is only small, especially under stable conditions. Emissions at such enclosed geometries may result in unrealistically high concentrations of primary compounds like NO x . With increasing grid resolution such narrow street canyons will be better represented on the 20 numerical grid, lowering the risk that high emissions coincide with regions where the flow is only poorly resolved. Hence, especially for lower grid resolutions, we recommend to carefully inspect the concentration output to treat high concentrations in the data analysis adequately.