Sensitivity of Pollutant Concentrations to the Turbulence Schemes of a Dispersion Modelling Chain over Complex Orography

: Atmospheric circulation over mountainous regions is more complex than over ﬂat terrain due to the interaction of ﬂows on various scales: synoptic-scale ﬂows, thermally-driven mesoscale winds and turbulent ﬂuxes. In order to faithfully reconstruct the circulation affecting the dispersion and deposition of pollutants in mountainous areas, meteorological models should have a sub-kilometer grid spacing, where turbulent motions are partially resolved and the parametrizations of the sub-grid scale ﬂuxes need to be evaluated. In this study, a modelling chain based on the Weather Research and Forecasting (WRF) model and the chemical transport model Flexible Air Quality Regional Model (FARM) is applied to estimate the pollutant concentrations at a 0.5 km horizontal resolution over the Aosta Valley, a mountainous region of the northwestern Alps. Two pollution episodes that occurred in this region are reconstructed: one summer episode dominated by thermally-driven winds, and one winter episode dominated by synoptic-scale ﬂows. Three WRF conﬁgurations with speciﬁc planetary boundary layer and surface layer schemes are tested, and the numerical results are compared with the surface measurements of meteorological variables at twenty-four stations. For each WRF conﬁguration, two different FARM runs are performed, with turbulence-related quantities provided by the SURface-atmosphere interFace PROcessor or directly by WRF. The chemical concentrations resulting from the different FARM runs are compared with the surface measurements of particulate matter of less than 10 µ m in diameter and nitrogen dioxide taken at ﬁve air quality stations. Furthermore, these results are compared with the outputs of the modelling chain employed routinely by the Aosta Valley Environmental Protection Agency, based on FARM driven by COSMO-I2 (COnsortium for Small-scale MOdelling) at 2.8 km horizontal grid spacing. The pollution events are underestimated by the modelling chain, but the bias between simulated and measured surface concentrations is reduced using the conﬁguration based on WRF turbulence parametrizations, which imply a reduced dispersion.


Introduction
Atmospheric transport processes over complex terrain are the result of a wide spectrum of flows acting and interacting at different spatial and temporal scales. Mountain chains can modify synoptic-scale flows, trigger mesoscale winds and alter the properties of turbulent fluxes. The inhomogeneity of the terrain and land cover modifies the mesoscale fluxes and influence the spatio-temporal evolution of microscale fluxes [1,2]. Recent applications of meteorological models over complex terrain show their ability to capture the main features during selected episodes, but the most local characteristics, such as slope winds, north-south in the first branch connected to the plain (Po Valley), and west-east oriented in the middle part, where Aosta town is located. The summer case study (25 August-1 September 2015) is representative of highpressure conditions (Figure 2a), triggering the development of a strong mesoscale circulation between the plain and the alpine valleys with higher intensity near the surface. The particulate matter concentrations follow the daily cycle of thermally driven winds, as shown by surface measurements and vertical profiles supplied by the ALC (Automated Lidar Ceilometer) in Saint-Christophe [17].

Meteorological and Chemical Stations
In the present analysis, twenty-four meteorological stations managed by the local bureau of meteorology (Centro Funzionale della Valle d'Aosta, https://cf.regione.vda.it/, last In the winter case study (24-31 January 2017), synoptic flows dominate transport processes, while weak winds and atmospheric stability prevail in the valleys (Figure 2b). The particulate concentrations rise concurrently with the flow rotation from the South-East, reaching surface values higher than 100 µg m −3 in Aosta, and do not decrease until the synoptic flow rotation from the North-West, as shown by the aerosol layer measured by the ALC [17].

Meteorological and Chemical Stations
In the present analysis, twenty-four meteorological stations managed by the local bureau of meteorology (Centro Funzionale della Valle d'Aosta, https://cf.regione.vda.it/, accessed on 22 December 2021) are chosen in the Aosta Valley to compare the meteorological simulations with observations. They are located at altitudes ranging from 341 m to 2500 m.a.s.l. (Figure 1). The considered quantities, at a time resolution of 30-min, are: temperature at 2 m (T2m), wind speed and direction at 10 m (W10m), relative humidity (RH), atmospheric pressure (P) and global radiation (SW). Hourly and daily time series of particulate matter of less than 10 µm of diameter (PM 10 ) and nitrogen dioxide (NO 2 ) concentrations come from five ARPA chemical stations, representative of different rural and urban areas. The complete list of the station characteristics is reported in Table 1. Table 1. List of the weather and chemical station used for comparisons. The measured quantities are shorted as follow: T2m-2 m-temperature, W10m-Wind speed and direction, RH-relative humidity, R-rain, P-atmospheric pressure, SW-global radiation, PM 10 -Particulate Matter less than 10 µ concentration, NO 2 -nitrogen dioxide concentration.  [19]. The middle d02 domain supplies the scale reduction necessary to simulate in detail the atmospheric fluxes in the innermost d03 domain that covers Aosta Valley and the neighboring areas ( Figure 3). The grid sizes and the features of each domain are summarized in Table 2. The WRF model nesting is applied in the two ways configuration. The vertical grid consists of 44 eta-levels, with 12 levels below 1000 m.a.g.l.

Station
In the innermost domain, the topography is the Shuttle Radar Topography Mission (SRTM, 3-arc-second resolution [20]), and the land use is the Corine Land Cover 2012 (CLC 2012, 100 m resolution, [21]) consisting of 44 classes. A detailed description of the model setup is available in [7].
Sub-grid processes are parametrized by applying specific schemes of the WRF environment. In all simulations, the longwave radiation scheme is the Rapid Radiative Transfer Model GCMs (RRTMG, [22]), the shortwave radiation one is the Dudhia [23] and the effects of topography shading on the shortwave radiation are considered over the d02 and d03 domains. The microphysics scheme is the WRF Single-Moment 6-class (WSM6, [24]), and the cumulus parametrization is activated for the d01 domain [25]. The Land Surface Model is Noah LSM [26]. The horizontal turbulent diffusion is parametrized with the Smagorinsky first-order closure, calculated in the physical space over d02 and d03 to correctly divide the horizontal and vertical fluxes over the slopes.
The three tested WRF configurations differ in the parametrizations of the vertical turbulent fluxes, and the PBL schemes will be used hereafter to indicate the WRF simulations ( Table 3). The set of PBL schemes are the Yonsei University (YSU, [27]) scheme, a diagnostic non-local scheme with an explicit description of the entrainment processes, the Mellor-Yamada-Nakanishi-Niino level 2.5 Eddy Diffusivity Mass-Flux (MYNN2.5 EDMF, [28]) scheme, with a mass-flux approach where the grid size limits the eddy size, making this scheme partially scale adaptive, and the Mellor-Yamada-Nakanishi-Niino level 3 (MYNN3, [29]) scheme, anturbulent kinetic energy prediction scheme with specific formulations of non-local turbulent transport. The formulation of the PBL height is based on the minimum heat flux level in the first-order scheme and the turbulent kinetic energy in the second-order schemes. The SL scheme set in the YSU configuration is the fifth-generation Mesoscale Model revised surface layer (MM5, [30]), while in the MYNN2.5 EDMF and MYNN3 configurations it is the Mellor-Yamada-Nakanishi-Niino surface layer (MYNN, [31]). In the MYNN2.5 EDMF configuration, the advection of the turbulent kinetic energy is activated over the d02 and d03 domains.  Microphysics WRF single-moment 6-class scheme [24] PBL physics + surface layer see Table 3 Land surface Unified Noah [26] Radiation RRTMG longwave scheme [24] Dudhia shortwave scheme [23] Topography  The applied air quality modelling system is based on the Flexible Air Quality Regional Model (FARM, http://farm-model.org, accessed on 22 December 2021), a 3D Eulerian model simulating dispersion and chemical reactions of atmospheric pollutants. The chemical transport model FARM is driven by an emission manager providing the emission rates over the domain and at the boundaries, a meteorological model describing the mean flow characteristics, and a micrometeorological pre-processor computing the meteorological information not available in the standard outputs of the applied meteorological driver, primarily the parametrizations of the sub-grid scale fluxes and the chemical-dependent parameters.
In the present work, the emission data over the FARM domain, included in the WRF d03 domain, come from the Aosta Valley emission inventory (managed by ARPA), and are pre-processed by EMMA (Emission Manager, http://doc.aria-net.it/EmissionManager, accessed on 22 December 2021) that prepares gridded hourly emissions. The initial and boundary conditions come from QualeAria (http://qualearia.it, accessed on 22 December 2021), a forecast system for the air quality in Italy and Europe. The meteorological driver is WRF, whose output fields (available every 30 min) are adapted by a module (GAP) to the FARM grid, featuring 32 vertical levels above the ground and 0.5 km horizontal grid spacing. The pre-processor SURFPRO (SURface-atmosphere interFace PROcessor, http://doc.aria-net.it/SURFPRO, accessed on 22 December 2021) is used to compute, based on the WRF wind, temperature, pressure, humidity, cloud cover and precipitation fields and land use information (CLC 2012): the turbulent and deposition fluxes of the chemical species described by FARM, through specific parametrizations [32]. In this work, the horizontal eddy diffusivity is calculated with the Smagorinsky first-order closure plus a stability-dependent term, and the vertical fluxes are described with a local approach where the diffusion coefficient depends on the friction velocity and the PBL height, which is set equal to a constant in the free atmosphere [33]. The PBL height formulation is based on the potential temperature profile in unstable conditions, while in stable conditions, it uses the friction velocity [34].
In this study, two modelling setups ( Figure 4) are tested for each WRF configuration (a total of six runs). In the first setup, the turbulence diffusion coefficients and scaling parameters are provided by WRF, allowing a better dynamical description since they are computed at each time-step taking account of the complete and updated set of meteorological quantities, while in the second setup, they are provided by SURFPRO driven by the 30-min WRF fields. The micrometeorological pre-processor provides the deposition velocities in both setups, as they are specific to the chemical species described by FARM.

Evaluation and Analysis
In order to assess the ability of the high-resolution meteorological model (in the three configurations) to reconstruct the atmospheric circulation during the episodes, the numerical results are compared with the surface measurements of the weather stations. The assessment of the model results is performed by looking at the time series of the individual stations to characterize the temporal evolution of the meteorological quantities for different altitude, slope and orientation conditions, and by grouping the data of all the stations by variable to derive a statistically robust evaluation of the model performance for the entire episode. The considered metrics are the Pearson's correlation coefficient (R 2 ), the Normalized Mean Bias (NMB), the Normalized Mean Gross Error (NMGE), the Root Mean Squared Error (RMSE) and the distribution quantiles, where: where M i is the i-th modelled value, and O i is the i-th observed value. Moreover, WRF temperature, water vapor mixing ratio and horizontal wind components are qualitatively examined with respect to the outputs of the COSMO-I2 model.
In order to assess the sensitivity of the simulated chemical concentrations to turbulence parametrizations, the FARM outputs are compared with the PM 10 and NO 2 observed concentrations. The comparison is performed between the modelled and the observed time series at the individual stations, to characterize the evolution of the pollutants in both rural and urban sites, and grouping the data of all the stations by variable to derive-if their number is sufficient-the same metrics as above. In order to interpret the different results of the two setups, the eddy diffusivities, the PBL heights and the friction velocities diagnosed by WRF and SURFPRO are compared. Furthermore, the PM 10 and NO 2 concentrations simulated by the WRF modelling chain are qualitatively compared with the outputs of FARM driven by COSMO-I2, to assess the impact of grid size in modelling the spatio-temporal evolution of concentrations over the considered area.  The performances of the applied WRF configurations over the entire episode are quantitatively similar. In particular, for temperature and relative humidity, the WRF YSU simulation features slightly higher correlations with the measures (Table 4), whereas for the wind speed, the WRF YSU, with the topography wind interaction activated, and the WRF MYNN2.5 performances are comparable. Indeed, the scale-adaptive mass-flux scheme for momentum activated in the WRF MYNN2.5 configuration contributes mainly in unstable summer conditions. The NMB, NMGE and RMSE confirm the agreement between simulated and observed surface temperatures and a general underestimation of the humidity ( Table 4).
The comparison between WRF and COSMO-I2 fields at 500 and 700 hPa shows a similar reconstruction of the synoptic circulation and of the mesoscale flows above the Alpine chain. At 850 hPa over the Po Valley, the models simulate a different evolution of the humidity, with higher values reconstructed by COSMO-I2, and of the thermally driven circulation from the plain to the mountains, with stronger up-valley winds simulated by WRF ( Figure 6).

Turbulence Parametrizations and Chemical Quantities
The comparison between the PBL heights, the friction velocities, and the eddy diffusivities computed by WRF and SURFPRO shows that, in stable conditions, the turbulence parameters diagnosed by both models are comparable, whereas, in unstable conditions, SURFPRO implies a stronger turbulent diffusion than the WRF configurations. Indeed, during the day, SURFPRO computes an higher PBL than the applied WRF configurations over the entire FARM domain. In Figure 7, the PBL height computed by the three WRF configurations (panels a-c) and by SURFPRO driven by the WRF YSU meteorological fields (panel d) are reported. Only results from SURFPRO driven by WRF YSU are shown, since the SURFPRO results derived with the other configurations of the meteorological driver are similar. The values and patterns are representative of the daytime conditions during the entire episode, characterized by a specific meteorological situation, and show that the PBL height spatial variability described by SURFPRO is reduced.  The eddy diffusivities are considered at five vertical levels ranging from 18 to 1792 m.a.g.l., and the major differences concern the horizontal eddy diffusivity near the surface, where SURFPRO diagnose values of one order of magnitude higher than the corresponding fields diagnosed by WRF, due to the contribution of the stability class-dependent term that dominates over the entire domain during the day. The vertical eddy diffusivity computed by SURFPRO features a spatial pattern similar to the friction velocity both near the surface and at elevated levels, while WRF provides a more pronounced vertical variability, with values of the order of magnitude of 10 1 m 2 s −1 near the surface and of 10 −2 m 2 s −1 at 480 and 1792 m.a.g.l. (one order of magnitude smaller than the corresponding SURFPRO field).
It is worth noting that the WRF YSU configuration reconstructs a strong gradient of the PBL height at the plain-mountain interface, highlighting that the schemes based on thermodynamic vertical profiles are not appropriate over complex terrain, where a multi-layer structure characterizes the PBL. Concerning the MYNN configurations, the PBL heights feature the same pattern of the turbulent kinetic energy: the shear energy production dominates in the valleys, whereas the buoyancy energy production is effective over the whole domain during the day (not shown).
The time series of the pollutant concentrations resulting from the two setups (see Figure 4) of the modelling system are compared with the measurements of NO 2 available in the air quality stations at hourly resolution, and of PM 10 available at the Aosta Piazza Plouves station at hourly resolution and at Donnas as daily averages. Both setups correctly reproduced the daily pollutant cycle, but the concentrations are underestimated during the entire episode. In the first setup, the bias between modelled and observed concentrations is smaller with respect to the second setup, due to the reduced dispersion implied by the WRF turbulence parametrizations. The NMB, NMGE, and RMSE derived grouping all NO 2 data confirm the better performance of the first setup (Table 5). The PM 10 hourly time series of the Aosta Piazza Plouves station confirms that the first setup reconstructs higher surface concentrations (Figure 8), and the daily averages of Donnas indicate a significant underestimation of the particulate matter near the South-East boundary of the domain (not shown). The latter result, and the presence of NO 2 and PM 10 concentrations remarkably higher than the seasonal averages in both rural and urban sites indicate that the transport of polluted air masses from the Po Basin, triggered by the plainmountain circulation, is an important contributor to the pollution episode. Moreover, the difference between observed and simulated concentrations, the latter following the daily cycle of the local emissions, can be explained by underestimated boundary emission rates. At a regional scale, the spatial patterns of the chemical concentrations resulting from the two setups of the WRF modelling chain are similar (Figure 9a-c). The comparison with the COSMO-I2 modelling chain (Figure 9d) shows that the spatial patterns reconstructed by the WRF modelling chain have a more significant variability both near the surface, where the flows are channelled in the tributary valleys), and in the free troposphere. In panels (a-c), the WRF modelling chain setup1 is depicted with the coordinate system is WGS84, while in panel (d), the COSMO-I2 modelling chain is depicted with coordinates in kilometers referred to the UTM (32N).

Meteorological Quantities
The modelled dynamic and thermodynamic variables' time series correctly reproduces the episode's atmospheric conditions. The relative humidity increases from 26 to 29 January, when the synoptic flow turns from West to South-East, and the wind speed in the stations located at high altitudes, representative of the synoptic circulation, is moderate during the same period ( Figure 10). On the contrary, both measured and observed wind speeds in the valleys are weak and irregular due to the atmospheric stability prevailing during the entire episode and triggering intermittent submeso flows. These irregular flows cause lower correlation coefficients between modelled and observed wind speed time series concerning the summer episode, characterized by a thermally driven circulation of moderate intensity near the surface. The thermodynamic variables feature higher correlations with measures, but a general underestimation of the humidity is present, with −0.17 < NMB < −0.14 ( Table 6). WRF and COSMO-I2 describe similar atmospheric circulations at 500 and 700 hPa, with flows blowing from South-East ( Figure 11) from 26 to 29 January and the advection of a warmer air mass at high altitudes on 27 January. At 850 hPa WRF describes in detail the impact of the Alpine chain, channeling and modifying the synoptic flows.

Turbulence Parametrizations and Chemical Quantities
The horizontal and vertical eddy diffusivities computed with SURFPRO are higher with respect to the corresponding coefficients computed by WRF at all considered vertical levels, ranging from 18 m to 1792 m.a.g.l., especially near the surface (Figure 12). In particular, within 100 m.a.g.l., the vertical eddy diffusivity diagnosed by the applied SURFPRO scheme is one order of magnitude higher than the corresponding WRF fields (on the order of 1 m 2 s −1 ) and features the same pattern of the friction velocity, which represents the major contributor to the turbulent dispersion during the episode due to the strong synoptic circulation and local stability. In contrast to the summer episode, the horizontal eddy diffusivities computed by the different models and schemes are of the same order of magnitude since the stability-dependent term is not dominant. The PBL height and the friction velocity, involved in the formulations of the diffusion coefficients, are higher on average with the SURFPRO formulation. The PBL height does not rise above 100 m.a.g.l. in the valleys during the entire episode due to the atmospheric stability and the consequent mixing inhibition, while the friction velocity is higher on the Alpine chain, where it reaches values close to 5 m s −1 , due to the shear of the synoptic flows. The spatial patterns of the TKE predicted by the MYNN configurations are similar to those of the friction velocity near the surface, where the shear energy production is dominant (not shown).  The time series resulting from the two setups are compared with the hourly data of NO 2 and PM 10 available at the air quality stations. In general, the daily cycle of pollutants is captured by the model, but the simulated concentrations do not rise significantly during the transport episode in both rural and urban sites ( Figure 13). The presence of high NO 2 and PM 10 concentrations in rural areas, and their persistence during the entire episode, may indicate a non-local component of the pollutant concentrations, mainly advected from the Po Basin by the synoptic circulation and trapped in the valley due to the strong atmospheric stability. The bias between modelled and observed concentrations is reduced with the first setup due to the reduced dispersion at local scale described by the applied WRF configurations. The NO 2 time series simulated by the first setup result in lower values of NMB, NMGE and RMSE (Table 7). Concerning the PM 10 concentrations (Table 8), the contribution of the Aosta Piazza Plouves station, where the PM 10 is overestimated, causes positive NMB and higher NMGE and RMSE in the first setup (Table 8) and highlights the need for long-term analysis to derive a robust statistical evaluation.
The concentrations simulated by the WRF modelling chain describe more realistically the impact of the orography channeling the advective fluxes from the Po Basin and the emission sources acting over the considered domain with respect to COSMO-I2 modelling chain (setup1 in Figure 14).

Conclusions
In this work, the meteorological and chemical quantities characterizing two pollution transport episodes, representative of typical atmospheric circulations over complex terrain, are reconstructed with a high resolution modelling chain based on WRF and FARM. Moreover, two modelling setups with different parametrizations of the sub-grid-scale fluxes and different coupling levels between the meteorological driver and the chemical transport model are tested.
The high-resolution WRF model, in the applied configurations, correctly reconstructs the meteorological circulations characterizing a summer and a winter episode, and triggering the transport of polluted air masses from the Po Basin to the North-Western Alps. Indeed, WRF surface fields have good correlations with the observations in both episodes. WRF and COSMO-I2 reconstruct similar synoptic circulations, and the higher resolution meteorological model describes in more detail the impact of the Alpine chain on the mesoscale circulation.
The first modelling chain setup, which uses the turbulent parameters directly evaluated during the WRF run, calculates higher surface concentrations of PM 10 and NO 2 . Indeed, the turbulent diffusion coefficients and scaling parameters resulting from the WRF configurations imply a reduced dispersion in both stable and unstable regimes with respect to the corresponding fields diagnosed by SURFPRO (second setup). The sensitivity to the modelling chain setups, featuring different formulations of the turbulence parametrizations, is greater than the sensitivity to the schemes set internally in WRF. The results of the first setup are in better agreement with the measurements and underline the added value of calculating the turbulence-related fields directly inside the meteorological model, instead of using external tools driven by a limited set of outputs of the meteorological model. The simulated pollutant concentrations are sensitive to the turbulence parametrizations at a local scale, where turbulent fluxes act, while the spatial resolution impact is primarily felt at the regional scale.
The difference between simulated and observed PM 10 and NO 2 concentrations, even if reduced in the first setup, is present using all the applied configurations of the modelling chain. It suggests an underestimation in the boundary and initial emission rates, particularly evident in these episodes characterized by a consistent advection from the South-East boundary of the domain (adjacent to the Po Basin). Purpose realistic boundaries and initial emission rates are fundamental to improving the modelling of such air pollution events.
In conclusion, the modelling of pollutant dispersion over complex terrain should take into account improvements in an accurate description of turbulence at the sub-kilometer scale, inside the "terra incognita", and in accurate initial and boundary conditions of the pollutant concentration fields.