Aerosol characteristics and polarimetric signatures for a deep convective storm over the northwestern part of Europe – modeling and observations

The Terrestrial Systems Modeling Platform (TSMP) was extended with a chemical transport model and polarimetric radar forward operator to enable detailed studies of aerosol–cloud–precipitation interactions. The model was used at kilometer-scale (convection-permitting) resolution to simulate a deep convective storm event over Germany which produced large hail, high precipitation, and severe damaging winds. The ensemble model simulation was, in general, able to capture the storm structure, its evolution, and the spatial pattern of accumulated precipitation. However, the model was found to underestimate regions of high accumulated precipitation (> 35 mm) and convective area fraction in the early period of the storm. While the model tends to simulate too high reflectivity in the downdraft region of the storm above the melting layer (mostly contributed by graupel), the model also simulates very weak polarimetric signatures in this region, when compared to the radar observations. The above findings remained almost unchanged when using a narrower cloud drop size distribution (CDSD) acknowledging the missing feedback between aerosol physical and chemical properties and CDSD shape parameters.

The kilometer-scale simulation showed that the strong updraft in the convective core produces aerosol-tower-like features, increasing the aerosol number concentrations and hence increasing the cloud droplet number concentration and reducing the mean cloud drop size. This could also be a source of discrepancy between the simulated polarimetric features like differential reflectivity (Z_DR) and specific differential-phase (K_DP) columns along the vicinity of the convective core compared to the X-band radar observations. However, the use of narrow CDSD did improve the simulation of Z_DR columns. Besides, the evaluation of simulated trace gases and aerosols was encouraging; however, a low bias was observed for aerosol optical depth (AOD), which could be partly linked to an underestimation of dust mass in the forcing data associated with a Saharan dust event.

This study illustrates the importance and the additional complexity associated with the inclusion of chemistry transport model when studying aerosol–cloud–precipitation interactions. But, along with polarimetric radar data for model evaluation, it allows us to identify and better constrain the traditional two-moment bulk cloud microphysical schemes used in the numerical weather prediction models for weather and climate.

1 Introduction

The effect of aerosol on clouds and precipitation through microphysical and radiative processes remains a major source of uncertainty in weather and climate prediction (Tao et al., 2012; Rosenfeld et al., 2014; Fan et al., 2016). In particular, improved understanding of the microphysical pathways of how aerosol affects cloud evolution (e.g., Rosenfeld et al., 2008; Koren et al., 2008; Yuan et al., 2011; Storer et al., 2014; Jiang et al., 2018; Igel and van den Heever, 2021) and precipitation (e.g., Rosenfeld, 2000; Stevens and Feingold, 2009; Shrestha and Barros, 2010; Li et al., 2011; Guo et al., 2016, 2018) is important for better prediction of extreme events. Many sensitivity studies using numerical models with various degrees of sophistication have been conducted to better understand these microphysical pathways with idealized/semi-idealized (e.g., Khain et al., 2005; Tao et al., 2007; Storer et al., 2010; Lebo and Seinfeld, 2011) or real data simulations (e.g., Noppel et al., 2010; Seifert et al., 2012; Morrison, 2012; Fan et al., 2013; Barros et al., 2018; Fan et al., 2018; Iguchi et al., 2020; Zhang et al., 2021; Trömel et al., 2021). While few of the above sensitivity studies have evaluated the model using radar reflectivity, polarimetric radar data which provide valuable information on cloud microphysical processes have not been fully exploited yet. In most of the numerical modeling studies, the aerosol physical and chemical properties have been held constant, and a large-scale perturbation of aerosol concentrations has been used for sensitivity studies. However, the classical assumptions made for “continental” or “marine” aerosols in the models do not reflect the actual local aerosol type, concentration, and its vertical profile or temporal evolution for any particular region on the globe. In fact, the meteorological settings, land cover, land use, and emissions strongly control the regional spectra of aerosol physical and chemical properties (e.g., Putaud et al., 2010; Martin et al., 2010; Shrestha et al., 2013). More recently, numerical modeling studies with a realistic aerosol distribution obtained by either downscaling region-specific aerosol profiles from a global aerosol model or using a meteorological model online coupled to a chemistry transport model have emerged (e.g., Rieger et al., 2014; Iguchi et al., 2020; Zhang et al., 2021). However, these studies have not fully exploited the potential of evaluating the model simulations against polarimetric radar observations. In this study, we use an online coupled meteorology–chemistry model (Baklanov et al., 2014), the Terrestrial Systems Modeling Platform (TSMP; Shrestha et al., 2014; Gasper et al., 2014) with the Aerosols and Reactive Trace gases (ART; Vogel et al., 2009) module for an ensemble simulation of a summertime deep convective storm over Germany at a kilometer-scale (convection-permitting) resolution. The main goal of the study is to (1) extend the TSMP with a chemistry transport model and a polarimetric radar forward operator to enable detailed studies of aerosol–cloud–precipitation interactions and their evaluation against polarimetric radar observations and (2) to demonstrate these new capabilities for a case study of a deep convective storm over Germany.

The paper is arranged as follows: Sect. 2 describes the observation data, model, forward operator (FO), and model setup used for the study. The first model evaluation of trace gases and aerosols with satellite and ground-based observations is presented in Sect. 3. The modeled aerosol physical and chemical characteristics during the storm event are presented in Sect. 4. The evaluations of modeled cloud microphysical processes and precipitation using polarimetric radar data are presented in Sect. 5. A detailed analysis of polarimetric features and aerosol characteristics is presented in Sect. 6. Finally, discussion and conclusions are presented in Sects. 7 and 8, respectively.

2 Data and methods

We study a summertime hail-bearing deep convective storm over northwestern Germany. The northeastward propagating storm was associated with the presence of pre-frontal convergence zones developed over this region on 5 July 2015. Scattered storms were prevalent throughout the day, with an isolated deep convective storm passing directly over Bonn from 15:00 to 16:00 UTC on 5 July 2015. Based on observations reported by the European Severe Weather Database (ESWD), large hail (2–5 cm in diameter) was observed over the Bonn region, including damaging lightning further north and heavy precipitation with severe wind (further northeast). A detailed discussion can also be found in Shrestha et al. (2022).

The study region encompasses the northwestern part of Germany bordering with the Netherlands, Luxembourg, Belgium, and France (Fig. 1). The region is characterized by multiple hills of the Rhine Massif, with heights ranging from 600 to 800 m and land cover including forest, agricultural land, and urban/rural area. The region also comprises extensive emissions by point (e.g., oil refineries and other industries) and area sources (e.g., extensive urban and rural areas, road transport, and extensive agriculture; Kulmala et al., 2011; Kuenen et al., 2014), making the region especially suited for this study. Additionally, due to the availability of the twin polarimetric X-band research radars in Bonn (BoXPol) and Jülich (JuXPol) and overlapping measurements from four polarimetric C-band radars of the German Weather Service (Deutscher Wetterdienst, DWD), along with the presence of the Jülich Observatory for Cloud Evolution (JOYCE; Löhnert et al., 2015), the region probably represents the best radar-monitored area in Germany.

Figure 1 Spatial pattern of topography over the study region. The Bonn radar domain is outlined in solid line, including the coverage of BoXPol and JuXPol (red circles). The dotted lines indicate the inner domain (excluding the relaxation zone) used to compute the domain average precipitation.

For coupled meteorology–chemistry modeling, we extend TSMP with the ART module. A forward operator is then used to transform the model outputs into radar space for evaluation with polarimetric radar observations from X-band radars. Available satellite observations and in situ observations are also used to evaluate the simulated trace gases and aerosols. A more detailed discussions about the observation data, model, forward operator, and the model setup are presented below.

2.1 Observations

In this study, we use polarimetric radar measurements of hydrometeors and ground/satellite-based estimates of aerosols and trace gases for model evaluation. The polarimetric radar measurements from the twin research X-band Doppler radars located in Bonn and Jülich (BoxPol and JuXPol; Diederich et al., 2015a, b) are used to investigate the microphysical characteristics of the deep convective storm. The polarimetric radar measurements provide valuable information about the horizontal reflectivity (Z_H), differential reflectivity (Z_DR), specific differential phase (K_DP), and cross-correlation coefficient (ρ_HV), which depend on hydrometeor shape, orientation, density, and phase composition and thus enable a detailed evaluation of the modeled microphysical and macrophysical processes. Z_H especially provides information on the size and, with that, on ongoing aggregation/riming processes. Z_DR mainly provides information on the shape of hydrometeors and does not depend on the number concentration, while K_DP is proportional to the concentration of hydrometeors. ρ_HV is mainly a measure of the hydrometeor diversity in the resolved radar resolution bin. Different patterns of these polarimetric variables enable us to identify ongoing cloud microphysical processes for precipitating systems. More comprehensive detail about radar polarimetry in general can be found in Ryzhkov and Zrnic (2019), Kumjian (2013), and Trömel et al. (2021), among many others. Further discussions about the polarimetric radar data are available in the Appendix (Sect. A1). Additionally, the RADOLAN (Radar Online Adjustment; Ramsauer et al., 2018; Kreklow et al., 2020) data from the German Weather Service (DWD, Deutscher Wetterdienst) is also used for model precipitation evaluation. RADOLAN is a gauge-adjusted precipitation product based on DWD's C-band weather radars available at hourly frequency in a spatial resolution of 1 km.

For the evaluation of modeled trace gases, remote sensing data from the Ozone Monitoring Instrument (OMI) aboard the Aura satellite is used. OMI provides valuable observations (e.g., O₃, NO₂, SO₂, and HCHO) to better understand the chemistry and dynamics of Earth's atmosphere. In this study, we make use of only OMI NO₂ v4.0 data (Krotkov et al., 2019; Lamsal et al., 2021a) to evaluate the spatial pattern of modeled NO₂ vertical tropospheric columns (VTCs). Other products were not used due to high uncertainty in their estimates for the timescales of evaluation used in this study. More discussions about the data are provided in the Appendix (Sect. A2).

Similarly, the Moderate Resolution Spectroradiometer (MODIS) 3 km aerosol product (MOD04_L3; Levy and Hsu, 2015) is used to evaluate the simulated spatial pattern of aerosol optical depth (AOD) at 550 nm. In addition, AOD level 2.0 (version 3) ground-based measurements from two AErosol RObotic NETwork (AERONET; Holben et al., 1998; Giles et al., 2019) stations over the domain are also used to evaluate the modeled AOD. These measurements have a better accuracy than MODIS (Giles et al., 2019) but are only available at a few locations.

2.2 Model

TSMP-ART v1.0 used in this study consists of the atmospheric model COSMO v5.1 (Consortium for Small-Scale Modeling; Steppeler et al., 2003; Baldauf et al., 2011) interfaced with ART v3.1 (Vogel et al., 2009), the land surface model CLM v3.5 (Community Land Model; Oleson et al., 2008), and 3D-distributed groundwater model ParFlow v3.1 (PARallel Flow; Ashby and Falgout, 1996; Jones and Woodward, 2001; Kollet and Maxwell, 2006; Maxwell, 2013). The three component models are coupled using the OASIS3-MCT coupler (Craig et al., 2017). COSMO-ART allows a comprehensive simulation of the two-way interaction between full gas-phase chemistry and aerosol dynamics with atmospheric processes (e.g., aerosol direct and indirect effects; washout of aerosols). Since ART v3.1 is already available as a module for the COSMO v5.1 model (which can be turned on with pre-processor flags), no extensive additional work was required to include the ART module in TSMP. As such, TSMP software was recently updated to include the ART v3.1 module with an extended version of the two-moment bulk microphysics scheme (Seifert and Beheng, 2006), including the hail class (Blahak, 2008, henceforth, SB2M). SB2M predicts the number and mass densities of cloud droplets, rain, cloud ice, snow, graupel, and hail, which are the zeroth and first moments of the particle mass distribution (PMD) that is assumed to follow a modified gamma distribution (MGD), as follows:

$$\begin{array}{}\text{(1)}& f\left(x\right)=N_{\mathrm{0}}{x}^{\mathit{\mu }}\exp (-\mathit{\lambda }{x}^{\mathit{\nu }}),\end{array}$$

with x being the particle mass, and parameters μ and ν determining the shape of the distribution. The specific hydrometeor mass q and specific number n can be derived by $q=Q/\mathit{\rho }$ and $n=N/\mathit{\rho }$, with ρ being the total density (air, vapor and hydrometeors).

The size–mass and velocity–mass relations of different hydrometeors are parameterized by power laws as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }D=a_g{x}^{b_g}\text{(3)}& {\displaystyle }{v}_T=a_v{x}^{b_v},\end{array}$$

with the (maximum) particle diameter D, terminal fall velocity v_T, and parameters a_g, b_g, a_v, and b_v.

The shape parameters μ and ν of the MGD remain constant for each hydrometeor class, and N₀ and λ can be diagnosed from the two prognostic moments. Furthermore, to mitigate the unphysical effects on the mean spectral particle mass $\stackrel{\mathrm{‾}}{x}=q/n$ coming from the separate advection and sedimentation of q and n, a minimum and maximum allowable mass limit is imposed for $\stackrel{\mathrm{‾}}{x}$ (x_min and x_max) at relevant places during the model time stepping. This is done by clipping n so that $\stackrel{\mathrm{‾}}{x}$ stays within [x_min,x_max]. For reference, all fixed parameters which were used in this study are summarized in Table 1.

Table 1 Parameters of the size–mass and velocity–mass relationships, following Eqs. (2) and (3) used in the SB2M. These refer to D in units of meters, x in kilograms, and v_T in meters per second. The last two columns contain the shape parameters of the assumed mass distribution. $D_{x,\min }=a_g{x}_{\min }^{b_g}$ and $D_{x,\max }=a_g{x}_{\max }^{b_g}$ are the diameters corresponding to the mass limits x_min and x_max and are added for better interpretation.

When the SB2M is coupled with the ART module, the cloud nucleation parameterization is based on the works of Fountoukis and Nenes (2005), Barahona and Nenes (2007), Kumar et al. (2009), and Barahona et al. (2010). Similarly, the ice nucleation parameterization is based on Barahona and Nenes (2009). A more detailed discussion about the implementation of the above nucleation parameterizations in ART is available from Bangert et al. (2012). The comprehensive activation parameterization works for a parcel of air containing an external mixture of soluble and insoluble aerosols. The activation rate is applied directly for newly formed clouds, while for existing clouds, the activation rate at the cloud base is calculated based on advection and the turbulent diffusion of particles into the cloud base (Bangert et al., 2012). Furthermore, for strong updrafts, in-cloud activation is also computed, for which the growth of existing cloud droplets is considered by assuming that they act as giant cloud condensation nuclei (CCN) that deplete supersaturation (Bangert et al., 2012). The activated aerosols as cloud droplets or nucleated into ice crystals are scavenged and removed. Besides the environmental and microphysical factors, the aerosol activation would also depend on its physical and chemical properties, which varies with elevation in the model. This can contribute to variable partitioning between the interstitial and activated aerosols as cloud droplets. Also, the parameterizations for the direct aerosol effect on the radiation and washout of aerosols by precipitation was turned on for the simulations. These formulations are all based on the prognostic aerosol population with 12 overlapping modes simulated in ART. Each mode is approximated by a lognormal distribution with uniform chemical composition across size. The 12 modes consist of the nucleation and accumulation mode for pure and mixed aerosol particles (sulfate, ammonium, nitrate, organic compounds, water, and soot), small, medium, and large particles for dust and sea salt, soot particles, and coarse particles (not used for the nucleation parameterization). These aerosol modes are coupled with gas-phase chemistry and strongly influenced by the atmospheric boundary layer evolution, advection, and anthropogenic emissions of gases and particles. An additional overview about the individual aerosol modes, chemical composition, and cloud interaction processes along with the aerosol dynamical processes can be found in Bangert et al. (2012).

For input of emission inventories, the online emission module developed earlier by Jähn et al. (2020) is used. This module makes use of pre-processed inventory data projected onto the model grid along with temporal and vertical scaling profiles for individual emission categories. A more detailed discussion about the pre-processing of emission inventories is presented in Sect. 2.4.

2.3 Forward operator

EMVORADO, the Efficient Modular VOlume scanning RADar Operator (Zeng et al., 2016) is COSMO's native radar forward operator. The FO uses model states and assumptions about the prescribed hydrometeor physical properties to compute the polarimetric radar variables, which are observed by X-band radars. FO requires consistency with the model, particularly regarding hydrometeor microphysics, i.e., size distributions and mass–size and velocity–size relations. For the online version run simultaneously with COSMO, this is ensured completely through variables shared between the modules. For an offline version run, this consistency is maintained manually. Here, we make use of the offline version, though, which is more flexible and allows us to re-run the FO with varied in-FO assumptions for, e.g., sensitivity analyses. More details about the FO is available in the Appendix (Sect. A3).

2.4 Model setup

The simulation is set up for an approximately 340 km × 340 km wide Bonn radar domain (Shrestha, 2021a) at a kilometer-scale resolution for the period 4 to 5 July 2015. For the initial and lateral boundary condition (IC/BC) of the atmospheric states in COSMO, data from the COSMO-DE ensemble prediction system (EPS; Gebhardt et al., 2011; Peralta et al., 2012) are used. The EPS data represent uncertainties in model physics and lateral boundary conditions by combining five model physics perturbations with four global models. An earlier study by Shrestha et al. (2022) showed that the statistics of the EPS are always stratified according to the four global models; i.e., the five members having the same global model are more similar to each other. In this study, we therefore only employ those five ensemble members that are based on the same global model of DWD (GME; Majewski et al., 2002). The ensemble simulation is used to reduce the uncertainty associated with meteorological forcings. The initial soil–vegetation states for CLM and ParFlow are obtained from spinups using offline hydrological model runs over the same domain (Shrestha, 2021a). For the initial and lateral boundary condition of trace gases and aerosols, we use the 6 h data from Model for Ozone and Related chemical Tracers, version 4 (MOZART-4; Emmons et al., 2010). The MOZART-4 data are available at a resolution of 1.9^∘ × 2.5^∘ with 56 levels. The COSMO model Processing Chain version 2.2 (available at https://github.com/C2SM/processing-chain, last access: 1 January 2022) was used for the pre-processing of the MOZART-4 data into ART variable states. This Python script maps the gases and aerosols (mass concentrations) from MOZART-4 to ART state variables. For the initialization of the number concentration of each aerosol mode, the default density and initial mode diameters in the ART module are used. Furthermore, we also assume that the aerosol has been in the atmosphere for a long time, where it could coagulate and mix, so 0.1 % and 99.9 % of the fine mode aerosols are assigned to mixed nucleation and accumulation mode, assuming a median diameter of these modes of 50 and 150 nm, respectively. The mapping from MOZART-4 to ART aerosol classes and the assumptions regarding median diameters are an additional source of uncertainty in the initialization of aerosols in the model.

The Copernicus Atmosphere Monitoring Service – Regional Inventory v4.2 (CAMS-REG v4.2; Kuenen et al., 2022) was used to prescribe the spatiotemporal emissions for the study. CAMS-REG v4.2 is a state-of-the-art gridded anthropogenic emission inventory developed for the European domain at a 0.1 ^∘ ×  0.05 ^∘ grid resolution, with a temporal coverage of 18 years (2000–2017). This emission inventory was pre-processed using the Python package “emiproc” (Jähn et al., 2020), available publicly through the C2SM GitHub organization (https://github.com/C2SM-RCM/cosmo-emission-processing, last access: 22 October 2022) for COSMO-ART variable states. First, the emission inventory data are projected onto the model grid, and then the temporal and vertical scaling profiles for individual emission categories are estimated. These inputs are then read during the model initialization, and the temporal and vertical emission profiles per category are applied online during the model run. In addition, the land cover data from Global Land Cover Map for 2009 (GlobCov 2009; Arino et al., 2012) is used for the biogenic volatile organic carbon (VOC) emissions. Furthermore, there is no emission of dust inside the model, and dust only comes from the MOZART-4 boundary conditions.

The ensemble simulation starts on 4 July 2015 at 06:00 UTC, and the model is integrated for 42 h. In all runs, a coupling frequency of 90 s is used between the atmospheric and hydrological components, which have a time step of 10 and 90 s, respectively. The model is integrated over the diurnal scale, and the output is generated at 5 min intervals and hourly intervals for evaluation with polarimetric radar data (only for a 3 h period) and aerosol measurements, respectively.

3 Evaluation of simulated trace gases and aerosols

First, the modeled trace gases and aerosols are evaluated with satellite- and ground-based observations. For comparison, the model data were also cloud screened. A threshold of 20 g m⁻² was used for the total column-integrated liquid and ice condensate for the cloud screening.

The vertical tropospheric column (VTC) is used to compare simulated NO₂ with satellite estimates from OMI. The VTC is an integral measure of the tracer from the surface to the tropopause. While it can be readily estimated from the model, the satellite estimates are dependent on the assumed vertical profiles of NO₂ in their algorithms. We acknowledge this uncertainty in the satellite estimates and the corresponding limitations of a direct comparison with the model data. However, it has to be stressed that this comparison is a very limited evaluation, since we are comparing the model with observations for a single day only. Both the satellite and the modeled VTC for NO₂ exhibit similar patterns, with relatively higher magnitudes over the northwestern lowlands and lower magnitudes over the Rhine Massif around 4 July 2015 at 12:00 UTC (Fig. 2). However, the model exhibits relatively higher magnitude of NO₂ over the foothills of the Rhine Massif near the emission sources (mostly from the mining regions and industry northwest of Bonn), which is not captured in the satellite retrievals. In order to compare simulated and observed NO₂ VTCs more quantitatively, we also interpolated the model output over the individual OMI pixels using an inverse distance squared algorithm (not shown here). The model exhibited a correlation coefficient of 0.46 with the observation.

Figure 2 Satellite (a) and model (b) estimates of integrated vertical tropospheric column (VTC) for NO₂ over the Bonn radar domain on 4 July 2015.

Figure 3 (a–b) Satellite and model estimates of aerosol optical depth (AOD at 550 nm) over the Bonn radar domain on 4 July 2015. The two available AERONET stations over the Bonn radar domain are also shown. (c–d) Time series of measured and simulated ensemble AOD over FJZ-JOYCE and MAINZ AERONET stations.

The modeled aerosol optical depth (AOD) is also compared with satellite retrievals from MODIS. In comparison with the MODIS data on 4 July 2015, the model tends to simulate relatively low AOD (0.1–0.3) over most of the domain (Fig. 3a–b). The MODIS data also show low AOD (0.1–0.3) over large parts of the domain but with pockets of high AOD scattered over the northern parts, which is not captured by the model. This bias in the modeled AOD can also be observed when comparing the modeled AOD with available AERONET station data over the region. The model generally tends to underestimate the AOD as estimated by the in situ measurements (when available). This is more prominent for the MAINZ station (Fig. 3d). However, within the spread of the ensemble members, the model also tends to the capture measured AOD over some period of times at Forschungszentrum Jülich (FZJ)-JOYCE station (Fig. 3c). In general, the above model evaluation with satellite data and ground-based measurements do build some confidence over the modeled aerosol and gaseous species.

4 Aerosol characterization

The modeled evolution of aerosol physical and chemical properties during the convective storm event are summarized in Figs. 4 and 5, respectively. Two different time periods are chosen as the storm propagates toward northwest, with strong updrafts in the synthetic sampling location at 15:00 UTC. The aerosol number concentrations of different modes (N_x) exhibit a strong variability in space and time. Figure 4a shows the spatial pattern of number concentrations of the sum of nucleation and accumulation mode for both pure and mixed aerosols (N_na) at 2 km height on 5 July 2015 at 14:00 UTC. At this time, the sampling location exhibits relatively low N_na compared to the western part, which has an extended patch of high N_na with an east to west extent. Over the next hour, this patch appears to be advected northwest, owing to the dominant southwesterly wind direction (Fig. 5a). At the same time, the spatial propagation of convective updrafts also plays a crucial role in lifting of aerosols to 2 km altitude. The evolution of the spatial pattern thus appears to be determined by a combination of horizontal advection and vertical updraft, with the latter additionally depending on the co-location with local emissions. Figure 4b shows the average aerosol size distribution for different modes and PM_2.5 (particulate matter with size <2.5 µm) chemical composition for a 9×9 grid cells box encompassing BoXPol at the center. At 2 km height, the dust particles dominate the aerosol mass, while soluble components make up only about 26 %. As expected, N_na is highest near the surface and generally decays with height. The magnitude of N_na is around 180 cm⁻³ at 2 km level. Also, a rightward shift in the aerosol size distribution of nucleation/accumulation mode can be observed that is associated with fresh aerosols near the surface and more aged aerosols in upper layers, with a larger mode at around 300 nm. The soot particles exhibit a multi-modal distribution with larger mode around 200 nm, while the dust particles exhibit a larger mode around 2000 nm. Figure 4c shows the meridional cross section of the aerosol number concentration (N_x) for a combination of different aerosol modes. As observed in the aerosol size distribution (Fig. 4b), N_na and N_soot exhibit higher concentration below 3 km height. Localized high values of N_na along the cross section are associated with local emissions. The dust aerosols exhibit a more horizontally homogeneous profile with a peak around 4 km, which is probably associated with a Saharan dust event. The multi-model forecast of dust from World Meteorological Organization (WMO) Barcelona Dust Regional Center (https://dust.aemet.es/products/daily-dust-products, last access: 22 October 2022) indicates the presence of Saharan dust for this particular event. The PM_2.5 concentration also shows peaks near the surface and near the melting layer but is associated with the N_na and N_dust, respectively.

Figure 4 (a) Spatial pattern of the aerosol number concentration for nucleation/accumulation (nuc./acc.; pure + mixed; N_na) at 2000 m height on 5 July 2015 at 14:00 UTC. The square with the “x” symbol at the center indicates the sampling location east of BoXPol, with the extent of a 9×9 box. (b) Average aerosol size distribution of different modes and PM_2.5 concentration for the 9×9 box. (c) Meridional cross section of the aerosol number concentration and PM_2.5 concentration passing through the sampling location. Also shown is the 0 ^∘C isotherm. (d) Ensemble vertical profile of aerosol number concentration for nuc./acc. (pure + mixed; N_na), soot (N_soot), dust (N_dust), and PM_2.5 concentration at “x”. (e) Ensemble time series of aerosol number concentration for nuc./acc. (pure + mixed) at 2000 m height at the BoXPol location. The blue and red lines correspond to times at 14:00 and 15:00 UTC, respectively.

Figure 5 Same as Fig. 4 but on 5 July 2015 at 15:00 UTC.

Figure 4d shows the area average vertical profile of the aerosol number concentration for different modes. The profiles are shown for the same 9×9 grid cells for five ensemble members. In general, all ensemble members exhibit similar profiles for this time period. Importantly, the aerosols exhibit a strong diurnal cycle owing to emissions, atmospheric boundary layer (ABL) evolution, updraft, and advection (Fig. 4e). At 2 km height, it generally peaks during the day and decays during the night (here only shown for nucleation/accumulation, nuc./acc., mode), except for periods with persistent convection or advection of aerosols. The situation at 15:00 UTC, when the aerosol distribution is strongly influenced by the deep convective event, is illustrated in Fig. 5. Due to the strong updraft associated with the convective storm, the aerosol size distribution of the nucleation/accumulation mode has become much broader (especially at 2–4 km height) as compared to the situation at 14:00 UTC. At the same time, the aerosol number concentration has increased, and the chemical composition (Fig. 5b) has changed significantly. The aerosol solubility and PM_2.5 mass concentrations at 2 km height have increased rapidly from 26 % to 46 % and 6.77 to 9.49 µg m⁻³, respectively. The simulated strong updraft over the sampling location also appears to generate localized aerosol towers reaching up to 15 km height (Fig. 5c). This increases the aerosol number concentration for all modes rapidly at higher altitudes (see Fig. 5d). We further discuss about these aerosol towers with polarimetric variables in following sections. In general, the variability in the location and magnitude of the simulated updraft associated with the convective storm produces the spread in the ensemble members.

5 Evaluation with polarimetric radar data

First, the modeled daily accumulated precipitation (5 July 2015) is evaluated with estimates from RADOLAN. Figure 6a–b show the spatial pattern of ensemble-averaged and RADOLAN-accumulated precipitation. In general, the model is able to capture the spatial pattern of the observed precipitation. However, the model underestimates the high precipitation in the northeastern part of the Bonn radar domain. This underestimation is also seen in the frequency distribution of the simulated and observed accumulated precipitation (Fig. 6c). While the domain average precipitation is similar to the RADOLAN data, all ensemble members tend to underestimate regions with high accumulated precipitation (> 35 mm). But, all ensemble members tend to slightly overestimate medium accumulated precipitation (10 to 30 mm).

Figure 6 Spatial pattern of accumulated precipitation. (a) Ensemble average from the model. (b) RADOLAN estimates. The black marker shows the location of BoXPol. (c) Frequency distribution of the simulated and observed accumulated precipitation. The inset shows the domain-averaged accumulated precipitation for each ensemble member (light gray color bar) and observation (black color bar) with 1 standard deviation (solid line above the bars).

Figure 7 Convective area fraction (CAF) of model ensemble members and observations. The two vertical bars define the time period used to compute contoured frequency altitude diagrams (CFADs) for the observation and model. The ensemble member, with a solid line, is used for polarimetric signature comparison. The CAF estimates from BoXPol or JuXPol are shown upon coverage and data availability. The gaps in the radar data represent times when the polarimetric signatures are strongly attenuated or if the storm extent is only partially covered by the radar.

The underestimation of high accumulated precipitation indicates that the model underestimates the high precipitation amounts associated with the core of the convective storm. This is also well seen in the time series of the convective area fraction (CAF; Fig. 7), which is estimated as the ratio of the storm area at 2 km above ground level (a.g.l. hereafter) with Z_H > 40 dBZ to the total area with Z_H > 0 dBZ. The masked storm area is generated using a storm tracking algorithm, which uses edge detection and overlapping areas between consecutive snapshots to track the storm. Observations from JuXPol and BoXPol exhibit high values of CAF in the early period of the storm (14:00 to 15:30 UTC), which is underestimated by all ensemble members. The ensemble members exhibit a similar pattern with increasing CAF after 15:30 UTC, when the simulated CAF matches the observed CAF more closely. However, such direct comparisons are always challenging due to mismatches in simulated and observed storm evolution in space and time, so we also conducted a qualitative exploratory analysis (using synthetic polarimetric variables at a lower altitude (1 km a.g.l.), near the melting level (4 km a.g.l.), and at a higher altitude (7 km a.g.l.) to find simulated convective storm structures closer in time and space to the radar observations. Based on this analysis, we compare the polarimetric signatures of the storm between one of the ensemble members (solid line; Fig. 7) and the BoXPol observations at 15:30 UTC.

Figure 8a) shows the plan position indicator (PPI) of polarimetric variables at 8.2^∘ elevation from BoXPol measurements. Near the melting level, the storm is characterized by high reflectivity (> 50 dBZ) and differential reflectivity (> 2 dB). At upper levels (beyond the convective core), the storm exhibits reflectivity in the range of 15 to 25 dBZ. The inflow region of the storm lies in the southeastern corner, which has relatively lower ρ_HV but high Z_H and Z_DR. The storm also exhibits an arc-like feature of high Z_DR along the eastern edge. Figure 8b shows the vertical cross section of the same polarimetric variables based on the gridded radar data along a north–south transect through the storm center. The convective core extends from −20 to 5 km, exhibiting high reflectivity (> 50 dBZ) from the surface up to 6 km height. A well-defined Z_DR column (> 2 dB) anchored to the surface and extending up to 6 km height is also visible along the cross section. Z_DR columns are distinct polarimetric signatures often found along the vicinity of the strong convective updraft core (Kumjian and Ryzhkov, 2008; Kumjian et al., 2014; Snyder et al., 2017b). K_DP columns (Ryzhkov and Zrnic, 2019; Snyder et al., 2017b) are also clearly distinguishable and co-located with Z_DR columns with slight inward offsets. High Z_DR and K_DP above the melting layer often indicate the presence of frozen raindrops, water-coated hail, and large size supercooled raindrops (Kumjian et al., 2014; Ryzhkov and Zrnic, 2019). Below the melting layer in the convective region, K_DP also has high magnitudes contributed by the melting of graupel/hail into raindrops. The high reflectivity values in the convective core with low ρ_HV (< 0.92) also indicates the dominance of the hail signature.

Figure 8 (a) Plan position indicator (PPI) plots of the horizontal reflectivity and differential reflectivity, specifically, the differential phase and cross-correlation coefficient at 8.2^∘ elevation measured by BoXPol on 5 July 2015 at 15:30 UTC. The dotted gray circles represent slant ranges for the chosen elevation angle, associated with heights of 1 km (lower level), 4.5 km (melting level), and 7 km (upper level). (b) Cross section of the same polarimetric variables from the gridded data. The vertical solid black line along the Y range in panel (a) indicates the location of the cross-sectioned plots.

Figure 9 (a, b) Model-simulated horizontal reflectivity and differential reflectivity, specifically, the differential phase and cross-correlation coefficient at the low-level (1000 m a.g.l.) and near-melting layer (4000 m a.g.l.) on 5 July 2015 at 15:30 UTC. The “x” symbol refers to the BoXPol location. The gray solid line indicates the location of the cross section. (c) Cross section of the same polarimetric variables. (d) Cross section of model-simulated hydrometeor density (QR is rain, QS is snow, QG is graupel, and QH is hail). Also shown is the 0 ^∘C line (solid black line) indicating the melting layer, contours of vertical velocity (5, 10, 20, 40 m s⁻¹) with QS, and contours of cloud ice density (QI) with QH.

Figure 9a–b show the spatial pattern of the synthetic polarimetric variables for the storm at 1 and 4 km a.g.l, as derived from the model simulation. Compared to the observations, the storm is already ahead of the BoXPol location but exhibits a similar structure compared to the observations. At lower levels (1 km a.g.l.), the storm exhibits an elongated zone with Z_H > 40 dBZ, which is also associated with relatively high Z_DR, K_DP but relatively lower ρ_HV. Near the melting level, the extent of the region with Z_H > 40 dBZ is much wider and also partly associated with high values of K_DP. However, relatively high values of Z_DR and lower values of ρ_HV are mostly constrained around the convective core. The meridional cross section of the synthetic polarimetric variables show that the storm top extends up to 13 km height, with an overshooting top up to 15 km height (Fig. 9c). The convective core also exhibits reflectivity > 50 dBZ up to 10 km height but is relatively narrow compared to the observation. A Z_DR-column-like feature protruding on top of the melting layer and anchored to the ground is also simulated; however, its magnitude is less than 0.8 dB. This is much weaker than the observed Z_DR columns with a magnitude of > 2 dB. Above the melting layer, Z_DR is generally weak (0 to 0.1 dB), with slightly higher values along the convective core. K_DP also exhibits relatively high values in the convective core, extending up to the storm top. ρ_HV is also relatively lower in the convective region and below the melting layer. In general, there is lack of polarimetric signal above the melting layer in the downdraft region of the storm, similar to an earlier study by Shrestha et al. (2022). The low variability and high values of synthetic ρ_HV can be attributed to the shortcomings in the FO assumption of the hydrometeor shape and orientation (Shrestha et al., 2022). The lack of a polarimetric signature in the downdraft region of the storm above the melting layer could be due to the deficiency in the FO for correctly modeling the scattering properties of snow and graupel which dominate this region, as discussed below.

The meridional cross section of the modeled hydrometeors shows the presence of supercooled raindrops in the strong updraft region, where the modeled vertical velocity above 8 km reaches 40 m s⁻¹ (Fig. 9d). The strong updraft also generates a warm anomaly above the melting layer (see the 0 ^∘C isotherm), below which rain is mainly produced by melting of graupel and hail. The melting of graupel and hail into raindrops produces the high K_DP below the melting layer. For ice hydrometeors, graupel dominates, with high-density surrounding the convective core. Graupel is responsible for the high reflectivity in the downdraft region of the storm above the melting layer. Cloud ice is located mostly above 8 km height and contributes to the high K_DP near the storm top. The self-collection of these ice particles leads to the formation of snow which extends further down to 6 km as it grows in size via aggregation. Hail mostly dominates in the strong updraft region of the storm with peaks in mass density adjacent to the supercooled raindrops. It also contributes to the high Z_DR values simulated in the convective region above the melting layer. The mean diameter of the supercooled raindrops is only around 0.1 to 0.3 mm, and the above-observed Z_DR-column-like signature is produced by the presence of water above freezing level due to melting of hail only. The mean hail size ranges from 0.1 to 13 mm (e.g., around 6 km height). During this time, the hail is also reaching the ground, starting from 15:25 to 15:40 UTC. In general, the Z_DR column usually appears 15–20 min before the hail reaches the ground (Ilotoviz et al., 2018). So, we also additionally explore this polarimetric feature in detail at earlier times in the following sections.

While the above analysis already indicates some shortcomings in the synthetic polarimetric signatures, the uncertainty due to mismatches between the space and timescales of synthetic and observed polarimetric variables also needs to be addressed by monitoring ensemble properties of the convective storm. So, additionally, the synthetic polarimetric variables from the ensemble simulations are compared to the observations from 14:45 to 15:30 UTC (see Fig. 7) using contoured frequency altitude diagrams (CFADs; Yuter and Houze, 1995) using the same extents and bin widths.

Figure 10a shows the CFADs of the polarimetric variables from BoXPol measurements. Z_H exhibits a narrow distribution at upper levels, with peaks around 15 to 25 dBZ. The distribution gradually broadens from the mid levels to the ground, with peaks around 15 to 40 dBZ near the melting layer. Z_DR has a narrow unimodal distribution above the melting layer with a peak around 0.14 dB. The distribution gradually broadens below the melting layer, with peaks shifting to 0.62 dB near the lower levels. K_DP also has a very narrow unimodal distribution with peak around 0.1 ^∘ km⁻¹. The distribution does exhibit a weak broadening from 8 km towards the surface. ρ_HV exhibits a broader distribution with peaks around 0.97 to 0.99 up to 10 km height. Above, the peak shifts towards 0.82 to 0.85.

Figure 10 CFADs of horizontal reflectivity and differential reflectivity, specifically, the differential phase and cross-correlation coefficient from 14:45 to 15:30 UTC on 5 July 2015. CFADs from the model are shown for five ensemble members.

Compared to the observations, the ensemble CFADs of synthetic Z_H exhibit a relatively broader distribution, with peaks around 20 to 35 dBZ in the upper levels. The peak of the distribution gradually shifts rightwards (30 to 40 dBZ) near the melting layer. Below the melting layer, the peak of the distribution shifts leftwards (3 to 25 dBZ), which also explains the lower CAF from the ensemble members compared to observations during this period. Similar to observations, synthetic Z_DR has a narrow distribution above the melting layer, with a peak around 0.11 dB. The distribution gradually broadens below the melting layer with additional peaks at around 1 and 2 dB. Similar to observations, synthetic K_DP has a very narrow unimodal distribution with a peak around 0.12 ^∘ km⁻¹. It also exhibits a weak broadening in the storm-top region and near the melting layer. For synthetic ρ_HV, the ensemble model CFADs show a very weak variability, with a peak at around 0.99 and a slight broadening below the melting layer.

6 Polarimetric feature and aerosol characteristics

The observations from X-band radar show the Z_DR and K_DP column as being one of the distinct features of this storm. However, the model is only able to simulate comparatively weak polarimetric features and contrasting aerosol-tower-like features. These polarimetric features and aerosol characteristics are therefore also explored for an earlier time at around 15:00 UTC. This time was chosen because it is also 25 min ahead of the hail reaching the ground and due to the availability of additional aerosol data (due to hourly output).

Figure 11a–b show the spatial pattern of Z_DR at 6 km height, along with the vertical wind speed. Enhanced Z_DR is present, surrounding a strong convective core. The width of the core is around 12 km, with the vertical speed exceeding 30 m s⁻¹. The forward flank downdraft and the rear flank downdraft is also visible. The meridional cross section shows the presence of a warm temperature perturbation above the melting layer in the convective core, which is mainly responsible for the melting of hail in the FO at relatively higher level, producing the ring-like Z_DR feature around the convective core.

Figure 11 Plan view (a) and vertical cross section (b) of aerosols, model states, and polarimetric variables. The plan views are shown at 6 km height, and all cross sections passing through the solid line are shown in the plan view. The 0 and −10 ^∘C isotherm is also shown in all cross sections. (a, b) Differential reflectivity (color fill) and vertical velocity (lines). The contoured solid/dashed red/blue lines indicate updraft and downdraft, respectively. The vertical wind speed contours are shown at the following intervals (−7.0, −5.0, −3.0, −1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0) in meters per second. (c, d) Rain and hail mixing ratios in filled and solid line contours, respectively. (e, f) Aerosol and cloud number concentrations in filled and dashed line contours, respectively. The cloud number concentration is contoured at interval of 500 cm⁻³, with a minimum of 100 cm⁻³. The aerosol concentration is shown for the pure and mixed nucleation and accumulation model aerosols.

The simulated Z_DR signal is mostly produced by the raindrops and hail particles (Fig. 11c, d). Raindrops dominate the convective core (above the melting layer) and downdraft region (below the melting layer) in terms of mass density. Hail mostly peaks northwest of the strongest convective core and also extends partly to the downdraft region above the melting layer. Above the melting layer in the convective core, the mean raindrop size is only around 0.1 to 0.3 mm, while the southwestern region exhibits grid-scale supercooled raindrops with size range of 1–3 mm, but part of its polarimetric signal is also masked by hail in the FO. The mean size of raindrops below the melting in the downdraft region is around 1 to 3 mm, which contributes strongly to the Z_DR signal besides the contribution from melting hail. The mean size of the hail particles is generally around 1 to 13 mm, with peak values around 6 to 9 km.

The comparatively small mean size of the hail particles and raindrops in the convective core could be due to the very high cloud drop number concentration simulated in the model (Fig. 11e, f). The cloud drop number concentration exhibits strong co-variability with the simulated nucleation/accumulation mode aerosol number concentration (N_na). The strong updraft increases the aerosol load in the convective core, which increases the aerosol number concentration and, consequently, the cloud drop number concentration, which varies from 100 to 3000 cm⁻³, leading to a very small size of cloud drops ranging from 5 to 25 µm.

7 Discussion

In this study, we extended the state-of-the-art terrestrial systems modeling platform with a chemistry transport module and a polarimetric forward operator. The model was then used to evaluate synthetic polarimetric signatures of a deep convective storm event over Germany with observations from X-band radar to better understand aerosol–cloud–precipitation interaction.

The model was also evaluated with satellite- and ground-based observations of trace gases and aerosols. The spatial pattern of NO₂ VTC was well captured by the model. This is consistent with an earlier evaluation of COSMO-ART by Knote et al. (2011), who also showed that the model was able to capture the spatial pattern and magnitudes compared to OMI estimates over Europe. Their study also showed that COSMO-ART underestimated the summertime AOD over much of Europe, compared to the estimates from MODIS, which is also consistent with the findings in this study. This indicates a possible model bias, which could be attributed to missing aerosol mass at lateral boundaries and inaccuracies in simulated aerosols within the domain (Knote et al., 2011). Additionally, the WMO Barcelona Dust Regional Center multi-model forecast shows dust AOD of 0.1 to 0.2 for this event; however, the model estimates of dust AOD are much lower, at around 0.04 to 0.06. This indicates that dust mass was possibly underestimated in the MOZART-4 data used in this study, which could be contributing the low bias of the simulated AOD.

In contrast to studies with fixed (e.g., climatological) aerosol distributions and properties, accounting for the full life cycle of aerosols using the ART v3.1 module introduces a strong diurnal cycle of aerosol physical and chemical properties which are modulated by synoptic winds and local convection. The typical large mode of the aerosol is around 300 nm, which is consistent with the assumptions made in SB2M runs (without the inclusion of ART v3.1 module). However, the number concentrations and chemical composition (hence solubility) of the aerosol exhibit strong variability in space and time. For example, during the convective event, the aerosol concentration and solubility at 2 km height rapidly increased. But, it has to be noted that the model could also be overestimating number concentrations near the ground, as found in an earlier model evaluation study by Knote et al. (2011) for many regions in Europe. The model simulation also shows a rapid increase in aerosol concentrations within the convective storm up to the overshooting cloud tops, generating aerosol towers with contrasting aerosol properties within and outside the storm. But, the uncertainty in the parameterization of the in-cloud processing of aerosols could also contribute to uncertainty in the simulated aerosol properties within the storm (Knote and Brunner, 2013).

In terms of accumulated precipitation, the model is able to capture the spatial pattern but underestimates the observed high precipitation amounts (> 35 mm) for all ensemble members. This finding is similar to results from an earlier study using TSMP with prescribed continental cloud nuclei (CN) and default ice nuclei (IN) concentrations (Shrestha et al., 2022). Also, similar to the finding in this study, the CAF is also underestimated in the early phase of the storm (14:45 to 15:30 UTC), compared to the observations. The underestimation of CAF could be associated with (1) reduction in collision and coalescence efficiency associated with small size of cloud droplets, (2) strong updrafts and high aerosol number concentrations, and (3) missing feedback between aerosol number concentrations and shape parameters governing the cloud drop size distribution. The kilometer-scale resolution of the current modeling study could be contributing to model-induced circulation enhancing the updraft speed (Poll et al., 2017, 2022), while the high aerosol number concentrations in the convective core resulting from the strong updraft contribute to a large number concentration of small cloud droplets.

In general, all ensemble members are able to capture the storm structure and evolution similar to the observations. However, the polarimetric signals above the melting layer are generally weak in the downdraft region, as also observed in earlier study (Shrestha et al., 2022), and also have higher reflectivity range compared to the observation. This is well captured in the CFADs compared to the observations. Above the melting layer, the model generally overestimates the horizontal reflectivity compared to the observation, which is primarily due to overproduction of graupel in the model. The predefined ice categories with fixed properties in a bulk microphysics scheme (e.g., SB2M used here) do not allow the simulation of full growth process for rimed particles like graupel or hail. This could contribute to the model bias in reflectivity in the downdraft region above the melting layer. A recent study by Milbrandt et al. (2021) has also shown that the three-moment representation of ice hydrometeors with the predicted particle properties (P3 scheme; Morrison and Milbrandt, 2015) improves the simulated reflectivity above the downdraft region for a hail-bearing storm.

Figure 12 (a) Modified gamma particle size distribution as a function of particle diameter (D_p) for the default and narrow cloud drop size distribution (CDSD). The bulk number concentration and mass density is 300 cm⁻³ and 1 gm⁻³, respectively. (b) Spatial pattern of ensemble-averaged accumulated precipitation for the default and narrow CDSD (solid contour lines with intervals at 20, 30 and 35 mm). (c) CFADs of horizontal reflectivity and differential reflectivity, specifically the differential phase and cross-correlation coefficient from 14:45 to 15:30 UTC on 5 July 2015. CFADs are shown for five ensemble members for narrow CDSD. The CFADs of the default experiment are shown in black contoured lines only.

Figure 13 Plan view (a) and vertical cross section (b) of differential reflectivity. The plan view is shown at 6 km height, and the cross section is passing through the solid line shown in the plan view. The 0 and −10 ^∘C isotherm are also shown in the cross section.

In terms of observed polarimetric features, the synthetic polarimetric variables also exhibit Z_DR-column-like feature (though of much weaker magnitude) along the updraft region, as in the observations. This difference may be attributed to a too small size of supercooled raindrops, but it may also be associated with the deficiency in the simulated updraft structure, recirculation of raindrops, and treatment of the slow freezing of raindrops (Kumjian et al., 2014; Snyder et al., 2015). Also, importantly, the Z_DR signal contribution from water-coated hail owing to wet growth process is missing. The current FO only has a parameterization for the melting of hydrometeors, but the water-coated hail particles due to wet growth are not included. Furthermore, the collision efficiency between frozen particle and supercooled droplets decreases with drop size, resulting in a weaker riming and hence producing smaller hail particles with lower fall velocity (Noppel et al., 2010). The study by Noppel et al. (2010), using the COSMO model with SB2M microphysics, showed that the continental CN concentration led to a weaker hail storm; however, an additional sensitivity study, by varying the shape parameters for cloud droplets producing narrow distribution, led to a different conclusion, indicating the missing feedback between the shape parameters of cloud droplets and CN concentrations. So, we also conducted an additional ensemble sensitivity study using a narrow cloud droplet size distribution (CDSD; see Fig. 12a). The parameters μ and ν, determining the shape of the distribution, were changed to 6 and 1, respectively, from the default value used in this study and referred to as narrow CDSD. With the narrow CDSD runs, all ensemble members still underestimate the CAF. The change in CDSD led to a delay in the onset of CAF evolution, with some ensemble members exhibiting relatively higher CAF. However, the CAF time series exhibits different variability for each ensemble member, suggesting the strong influence of lateral boundary conditions. The domain average precipitation and the spread of the frequency distribution of the precipitation are similar to the default runs. Only the spatial location of high precipitation for the ensemble average in the northeastern part of the domain is slightly shifted (Fig. 12b). However, the narrow CDSD does show improvement in the simulated Z_DR-column-like features, which is more well defined than the default experiment with larger mean raindrop size (0.5–1 mm) above the melting layer (Fig. 13). Besides, the narrow CDSD does affect the CFADs of the storm in terms of polarimetric variables (Fig. 12c). At upper levels, the peaks of Z_H shift to higher magnitudes at 20 to 35 dbZ. And, the distribution gradually broadens, and the peak shifts rightward (30 to 40 dBZ) near the melting layer. Below the melting layer, the distribution shifts rightward, as simulated before (with default CDSD), with peaks around 5 to 25 dBZ. CFADs of Z_DR also exhibit multimodal distribution below the melting layer with additional peaks around 0.87 and 1.87 dB, which is slightly lower than the default CDSD runs. Above the melting layer, the peak of Z_DR remains at 0.11 dB. The CFADs of the K_DP and ρ_HV also exhibit a similar peak around 0.11 ^∘ km⁻¹ and 0.99, respectively. But, in general, there is an increase in the spread of all the polarimetric variables. This could probably indicate the importance of the shape parameters of the hydrometeors for improving the simulated polarimetric signature of the storm.

8 Conclusions

While acknowledging the model biases and uncertainty in the simulated aerosol properties, the inclusion of prognostic aerosol is a way forward for a better understanding of the aerosol–cloud–precipitation interactions. During the convective storm event, the model generates aerosol-tower-like features with contrasting physical and chemical properties compared to the background.

At diurnal scales, the model is able to capture the spatial pattern of the precipitation; however, the comparison with the polarimetric observations indicates possible deviation in the ice hydrometeor partitioning above the melting layer (especially in the downdraft region of the storm), the size of supercooled raindrops and hail in the vicinity of the convective core, and the mechanism of rain production below the melting layer – hence the particle shapes and concentration. Constraining the CDSD did not produce any drastic difference in the simulated precipitation but produced an improvement in the simulated Z_DR-column-like features. Thus, running simulations with prognostic aerosol and the use of a forward operator can also help to constrain the cloud droplet size distribution in the model. Besides the shortcomings in the traditional two-moment bulk scheme used in this study, the effect of the model grid resolution and its impact on the structure of the storm updraft and the effect of a simulated high number of aerosol concentration which is lifted in the convective core, and hence the polarimetric signature in the vicinity of the convective core, can also not be neglected.

Thus, future aerosol–cloud–precipitation interaction studies using models should make an effort to include prognostic aerosol models and evaluate the cloud microphysical processes using polarimetric radar data to identify and improve the cloud microphysical parameterization in the current numerical weather prediction model used for weather and climate prediction.

Appendix A

A1 Polarimetric radar data

The X-band Doppler radars are operating at a frequency of 9.3 GHz, with a radial resolution of 100–150 m and a scan period of 5 min. Both radars produce volume scans at different elevations, mostly between 0.5 and 30^∘. These volume scans are also used to interpolate the polarimetric radar data from the native polar coordinates to Cartesian coordinates at 500 m horizontal and vertical resolution, using a Cressman analysis with a radius of influence of 2 km in the horizontal and 1 km in the vertical. A threshold of 0.8 in ρ_HV was imposed on the gridded data to ensure that clutter is filtered out without removing useful meteorological information.

The polarimetric variables Z_H and Z_DR are potentially affected by radar miscalibration, partial beam blockage, and (differential) attenuation, especially at smaller wavelengths (C band and X band), and their correction especially in deep convective, hail-bearing cells gives rise to additional uncertainties (e.g., Snyder et al., 2010; Shrestha et al., 2022). Although K_DP estimates are not affected by miscalibration and attenuation, they can be substantially affected by the uncertainty in the quantification of the backscatter differential phase (δ), which is particularly important when hydrometeor sizes are in the range of, or larger than, the radar wavelength (Trömel et al., 2013). A more detailed discussion about the calibration, clutter filtering, and attenuation correction of the polarimetric radar data can be found in Shrestha et al. (2022). It is important to note that errors in the estimates of polarimetric radar variables might arise due to the assumptions made in the attenuation correction algorithm and due to uncertainties in the contribution of the backscatter differential phase to the total differential-phase shift. We acknowledge these limitations in the study and concentrate more on patterns and not so much on the actual magnitudes of the polarimetric moments.

A2 Trace gases and aerosols

NO₂ is a key anthropogenic air pollutant and precursor of aerosols. The product comes with an estimated uncertainty for each VTC that could be used for the comparison between satellite and model. According to Boersma et al. (2011), Boersma et al. (2018), and Lamsal et al. (2021b), typical uncertainties are of the order of 30 % under clear-sky conditions. This should also hold for the data used in our study. The OMI estimates of VTC NO₂ are filtered for data points with VcdfQualityFlags = 0 and CloudRadianceFraction < 0.5 (clear sky data). We did not include OMI HCHO and O₃ for the following reasons: the HCHO product is extremely noisy. The retrieval uncertainty is 50 %–105 %, with the lower end being valid only for highly polluted locations. HCHO products are therefore usually only presented as monthly, seasonal, or yearly averages (e.g., De Smedt et al., 2021). The same argument holds for SO₂. Comparing O₃ would be quite interesting, but there is no official OMI tropospheric O₃ product. Also, comparing the total column O₃ from OMI would not be meaningful, as the column is strongly dominated by the stratosphere. Note also that, due to the long lifetime of O₃, we would expect only very small gradients in the model domain.

A3 Polarimetric radar forward operator

For an offline run, EMVORADO requires as input the atmospheric fields of mass and (for SB2M) number concentrations of the six hydrometeor classes (cloud liquid, rain, cloud ice, snow, graupel, and hail) of temperature and of the three wind components. Other parameters that affect forward-modeled polarimetric radar observables are insufficiently constrained by the COSMO model, and assumptions need to be made within the FO. This regards, e.g., the phase partitioning of hydrometeors during melting, the shape and orientation of particles, and the heterogeneous microstructure of frozen hydrometeors.

Like essentially all bulk scheme models, SB2M does not provide a prognostic melt fraction, and hydrometeors are either (completely) frozen or liquid. All meltwater is assumed to be shed instantaneously and transferred into the rain hydrometeor class; hence, no mixed-phase hydrometeors are predicted. Liquid water and ice exhibit significantly different dielectric properties in the radar frequency region, which leads to strong changes in the reflectivities where a phase change takes place. The melting layer is hence appearing very prominently in radar observations, particularly in stratiform situations, as layer of enhanced reflectivity known as the radar bright band. In order to be able to simulate such features, the forward operator needs to employ a melting scheme that predicts the occurrence of mixed-phase, wet frozen hydrometeors based on the single-phase model hydrometeors. EMVORADO employs a melting scheme that assumes a certain fraction of the frozen hydrometeor mass to be liquid (in contrast to, e.g., Wolfensberger and Berne, 2018, and Jung et al., 2008, who redistribute a part of the rainwater back into the frozen hydrometeor classes, i.e., unshed some rainwater). EMVORADO models the liquid water fraction dependent on the size of the hydrometeors (considering that small particles melt faster than large ones) and the ambient temperature T (Blahak, 2016). Wet hydrometeors start to occur when T exceeds a threshold T_meltbegin and are assumed to be completely melted when T_max is reached, where T_max by default is determined dynamically from the model hydrometeor field and T in the vertical column. Setting T_meltbegin accordingly allows for wet frozen hydrometeors at sub-zero temperatures, covering the case of the upward transport of liquid water and wet hydrometeors that do not (re-)freeze instantaneously in convective updrafts. Through the temperature dependence parameters, which are specific to each frozen hydrometeor class, the melting scheme can be adjusted by the user. Unless noted otherwise, in this study we apply EMVORADO's default melting scheme parameters (see Table A1).

Table A1 Overview of EMVORADO melting scheme setup used in this study.

The shape and orientation of the hydrometeors significantly affect the polarimetric radar parameters but are entirely unconstrained by the COSMO model. Here we make use of the polarimetric mode of EMVORADO, which so far applies the T-matrix scattering method for one- (Mishchenko, 2000) or two-layered (Ryzhkov et al., 2011) spheroidal particles. All hydrometeors are assumed to be oblate spheroids (except for liquid cloud particles that are modeled as spheres using Rayleigh scattering) with hydrometeor class-specific and size- and melt-fraction-dependent parameterizations of shape and orientation of the hydrometeors, as given by Ryzhkov et al. (2011). The effect of orientation distributions is considered using the angular moments approach outlined in Ryzhkov (2001); Ryzhkov et al. (2011, 2013).

In order to allow the fast calculation of the radar observables, lookup tables of bulk scattering properties are pre-calculated, tabulating basic (additive) quantities per hydrometeor class over bulk (mean) mass, temperature, and melting T_max. These are then added up over the six hydrometeor classes and converted into the polarimetric radar observables. Beside the reflectivity factor in horizontal polarization Z_H, in this study we focus on the differential reflectivity Z_DR, the co-polar cross-correlation coefficient ρ_HV, and the specific differential phase (K_DP). In short, Z_DR is the difference in the (log or dBZ space) reflectivities in the horizontal and vertical polarization, ρ_HV, the correlation between reflectivities in horizontal and vertical polarization within the measurement volume, and K_DP the phase difference between the horizontal and vertical polarized wave returns. A more comprehensive description can be found, e.g., in Kumjian (2013).

EMVORADO is capable of simulating the sensing process, including scanning, beam tracing, beam blockage, beam pattern, attenuation. This allows us to directly simulate observation equivalents like 3D volume scans. However, here we make use of the radar parameters calculated on the model grid, i.e., neglecting any sensing effects. Further details about the FO and its sensitivity to assumed parameters for the hydrometeors can be found in Shrestha et al. (2022),Trömel et al. (2021), and Shrestha et al. (2022).