High Resolution Chemistry Transport Modeling with the On-Line CHIMERE-WRF Model over the French Alps—Analysis of a Feedback of Surface Particulate Matter Concentrations on Mountainous Meteorology

: Air pollution is of major concern throughout the world and the use of modeling tools to analyze and forecast the pollutant concentrations in complex orographic areas remains challenging. This work proposes an exhaustive framework to analyze the ability of models to simulate the air quality over the French Alps up to 1.2 km resolution over Grenoble and the Arve Valley. The on-line coupled suite of models CHIMERE-WRF is used in its recent version to analyze a 1 month episode in November–December 2013. As expected, an improved resolution increases the concentrations close to the emission areas and reduced the negative bias for Particulate Matter that is the usual weakness of air quality models. However, the nitrate concentrations seem overestimated with at the same time an overestimation of surface temperature in the morning by WRF. Different WRF settings found in the literature are tested to improve the results, particularly the ability of the meteorological model to simulate the strong thermal inversions in the morning. Wood burning is one of the main contributor of air pollution during the period ranging from 80 to 90% of the Organic Matter. The activation of the on-line coupling has a moderate impact on the background concentrations but surprisingly a change of Particulate Matter (PM) concentrations in the valley will affect more the meteorology nearby high altitude areas than in the valley. This phenomenon is the result of a chain of processes involving the radiative effects and the water vapor column gradients in complex orographic areas. At last, the model conﬁrms that the surrounding glaciers are largely impacted by long range transport of desert dust. However, in wintertime some outbreaks of anthropogenic pollution from the valley when the synoptic situation changes can be advected up to the nearby high altitude areas, then deposited.


Introduction
The deterioration of air quality is one of the major environmental threats throughout the world. The impact of pollutants on health, ecosystems, climate and renewable energy production are clearly highlighted by many studies [1][2][3][4][5][6][7]. Regarding health, statistically significant positive correlation between the number of influenza-like illness cases and Air Quality Index were found in a recent study in China [8]. Often called short-lived climate pollutants (SLCPs), black carbon-a component of Particulate Matter (PM)-tropospheric ozone and methane contribute to both the warming of the climate as well as air pollution.
Modeling tools, the chemical transport models (CTMs) such as CHIMERE, are useful to evaluate, analyze and forecast the air quality from urban to global scales [9]. In France, the national forecast platform PREV'AIR [10] provides a 4-day air quality forecast now at approximately 4 km resolution over France and at approximately 10 km over Europe. High pressure systems developing over Western Europe from the north of France to Scandinavia lead to favorable conditions for the development of massive PM pollution events over Western Europe. These kinds of conditions are frequent from fall to spring and have been described in previous studies [11][12][13]. A large amount of ammonium nitrate is produced by the reaction of ammonia (mainly emitted by the agricultural sector) and nitric acid issued from the oxidation of nitrogen oxides (mainly associated to road transport and industrial activities).
The mesoscale meteorological model WRF [14] is widely used to simulate local meteorology and shows strong limitations over hilly areas for instance in the case of Californian valleys [15,16] during highly stable periods. Some case studies are reported in the literature investigating the complex phenomena of dispersion and transport of pollutants involving valley and slope winds [17]. The WRF-chem model was used over the Alps in Switzerland at 2 × 2 km resolution [18], and they showed the ability of the model to capture the temporal evolution of three-dimensional data for a variety of air pollutants and meteorological parameters, but a post-processing statistical model was developed to improve the raw model outputs. Results with WRF-chem in [19] also clearly highlight the difficulties of capturing meteorological parameters in the complex terrain of the Kathmandu Valley and reproducing subgrid-scale processes with a horizontal resolution of 3 × 3 km 2 .
In [20], the effect of wind dynamics on the dispersion of passive scalars is identified by tracking areas prone to stagnation, recirculation and ventilation. In a companion paper [21] they show the impact of a valley-wind system on the transport of passive tracers in the stably-stratified atmosphere of a valley dynamically decoupled from the atmosphere above. The amplitude of the oscillations can be as high as 50% of their mean value, implying that averaged values in an urbanized valley may disguise high instantaneous and potentially harmful concentration values. Our work is also an opportunity to evaluate at local scale the impact of PM emissions in the valleys on the meteorological conditions of nearby clean air altitude areas through the modification of radiative budget. So far this topic has been often addressed at regional or global scales with a focus on remote areas like the Artic [22,23].
In this work the CHIMERE model in its most recent version on-line coupled with WRF [24,25] merged with developments on the organic chemistry as presented in [26] is used to simulate a long wintertime period over the Alps covering several PM episodes. To our knowledge this study is a first published work dealing with air quality simulation over French Alpine valleys at resolutions close to 1 km. The selected period is November-December 2013 where several episodes was identified over our two target domains: Grenoble and the Arve Valley. The scope of the paper is fivefold: • Analyze the transport of pollutants under winter conditions in Alpine valleys up to the glaciers through the wet and dry deposition processes.

Model Set-Up
The CHIMERE configuration is summarized here, but the reader can refer to the reference CHIMERE publications [9,[24][25][26] for details on the corresponding model components and references as well as non user-specific model characteristics. The on-line coupling was developed between CHIMERE and WRF following two phases: the direct effect by [24] and the indirect effect in [25]. The aerosol direct and indirect effects were studied over west Africa in the summer of 2016 using the coupled WRF-CHIMERE regional model including aerosol-cloud interactions [27]. The gas-phase chemical mechanism is MELCHIOR2, which consists of a simplified version (more than 40 species and 120 reactions) of the full chemical mechanism based on the concept of chemical operators. Modeled particulate matter includes primary particulate matter, and secondary inorganic (nitrate, sulfate, ammonium based on the ISORROPIA thermodynamic equilibrium model) and organic aerosol resulting from the oxidation of the relevant anthropogenic and biogenic precursors and gas-particle partitioning of the condensable oxidation products [26]. Biogenic emissions are computed with MEGAN version 2.1 [28], sea-salt and mineral dust emissions from desert and agricultural areas are also considered.
The particle size ranges from 10 nm to 40 µm over 10 bins. In this paper SIA, SOA, OM, EC, DUST, SALT, PPM respectively refer to Secondary Inorganic Aerosol (sum of nitrate, sulfate and ammonium), Secondary Organic Aerosols (anthropogenic and biogenic in origins), Organic Matter, Elemental Carbon, natural mineral dust, sea salt and Primary Particle Matter (primary anthropogenic carbonaceous and non-carbonaceous species). Horizontal transport is solved with the second-order Van Leer scheme [29]. Subgrid scale convective fluxes are considered. Once the depth of the boundary layer is computed, vertical turbulent mixing can be applied by following the K diffusion framework, which follows the parameterization of [30], but without counter-gradient term. A minimum vertical eddy diffusion (KZZ) is set to 0.1 m 2 s −1 and a minimum Planetary Boundary Layer (PBL) height of 20 m is applied to avoid unrealistic concentrations.
For the meteorology, the WRF 3.7.1 version is used. The WRF simulation has been nudged to National Centers for Environmental Prediction (NCEP) final analysis meteorological fields Global Forecast System (GFS) at 1 × 1 • . and 6-hourly resolution at the coarsest model initial and domain boundaries (ds083.2 dataset, [31]).
For the reference simulation performed in this study, REF01, the main WRF settings are summarized in Table 1 [32][33][34][35][36][37]. The choice for the Planetary Boundary Layer (PBL) scheme YSU is supported by [38] and proposed as a "best practice" by the WRF development team as well as in [39]. Overall, the use of YSU exhibits good performances but as shown by [40] there is no universal good schemes, and a part of this study will be to test other configurations in this peculiar region. For high resolution simulations the common practice is to avoid to use a parameterization of the Cumulus cu_physics=0. The PBL computed by WRF is directly used in CHIMERE with a minimum value of 20 m to avoid unrealistic overshoot on the computation of concentrations.
For the anthropogenic emissions, the main pollutant's emission values (PPM2.5, PPMcoarse, NOx, NH3, CO, NMVOC, SO2) are from EMEP (Co-operative program for monitoring and evaluation of long range transmission of air pollutants in Europe) program [41]. For France, the emission inventory benefits from a spatial improvement using 1 × 1 km proxies for the main activity sectors. Similarly as described in [9], the proxy is directly calculated as the value of the yearly emissions for a given pollutant and activity sector, then this proxy is used to reallocate the emissions from the national level to the final target grid. These accurate data are from the French National Inventory accessible on request [42]. Organic Matter is split into several Primary Organic Aerosol (POA) species partitioned between the gas and aerosol phases which react with oxidant species. The scheme is described in [26] as well as the split and scaled factors to account for Intermediate and Semi-Volatile Organic Compunds (IVOC/SVOC) emissions and gap filling on activities particularly on residential wood burning emissions as described in [43].
At the boundaries of the European domain, monthly climatologies from global model simulations are used. In this study, outputs from LMDzINCA [44] are used for all gaseous and aerosols species, except for mineral dust where the simulations from the GOCART model are used [45]. In the reference simulation the model is on-line coupled without feedback of the chemistry on the meteorology.

Domains
Grenoble is a medium city located in the Alps at around 200 m altitude. The city is surrounded by massif peaking at 2000 to 3000 m on the East (Belledone, Taillefer), North (Chartreuse) and South-West (Vercors). The Arve valley is a mountain valley in the French Alps, ranging from the Col de Balme (2191 m, non-paved pass on the French-Swiss border) to the surroundings of Bonneville, downstream. This valley is between Geneva and the Mont-Blanc tunnel (close to Chamonix), a major itinerary for both light and heavy-duty vehicles between France and Italy, as well as a secondary road connecting France with Switzerland at the Col des Montets, paved pass at 1461 m a.s.l. The valley is surrounded on both sides by mountains higher than 3000 m, including the Mont Blanc Massif to the South-East. (Figure 1). These two target areas "Grenoble" and "Arve Valley" suffer from frequent pollution episodes as depicted in the introduction. For this study, four different domains are defined and explained in Table 2: EUR01 over Western Europe, ALP0033 for the whole French Alps, GRE0011 and ARV0011 to have fine zooms over Grenoble and Arve Valley. The target domains GRE0011 and ARV0011 have a horizontal resolution close to 1.2 km, they are nested in a one-way strategy from the continental scale (EUR01) as depicted in Table 2 with an intermediary domain ALP0033. All domains are regular in longitude and latitude.  While a spectral nudging of the horizontal wind components, temperature and the specific humidity is applied to avoid any divergence for the coarser domain, this option is turned off for the nested domains to fully benefit from the WRF physics at local scale. In the CHIMERE-WRF suite the 2-way nesting is not possible then for the 1-way nesting the ndown.exe WRF utility is used to provide the necessary inputs for the child domain from the parent domain as reported in recent studies [47,48] and the WRF documentation [49].
For all domains, the vertical structure of eta levels for WRF is the same with 33 levels (eta_levels) until 50 hPa as: 1 [36] as described by [50]. The first CHIMERE pressure coordinate is 998 hPa, the top pressure is 300 hPa. The first three level top heights at the highest resolution are respectively about 15, 35 and 60 m.

Period
Several pollution episodes occurred in France in November and December 2013 affecting the whole country and particularly over the Northern and Eastern France. Over the Alps, the weather conditions were often very stable with frequent High Pressure systems developing from North of France to Central Europe. As usual in these situations, Low Pressure systems over the north of Italy and the Mediterranean basin can develop and some precipitations were recorded on short periods over the Alps, for instance in the Grenoble area: Some synoptic maps of the meteorological situation are provided in Supplementary Material Figure S1 for these key dates. The official air quality forecasting system PREV'AIR [10] archives the observed and modeled dataset and is accessible at www.prevair.org. Several episodes were recorded from 28 Novermber to 19 December, this work focusing on the 27-30 Novermber 2013 (as EPI) and 11-18 December 2013 (as EPII) episodes. Under these synoptic meteorological conditions during this period, no specific desert dust outbreaks from the Sahara was recorded, a background of dust concentrations affects the regions. During this season in the Alps, no local emitted dust was simulated since the soil humidity is too high associated with too low wind speed. Therefore, a period from 10 November to 20 December was simulated but only the period from 15 November to 20 December is analyzed, the first 5-day block is used for as a spin-up for initialization. In many locations, the EU daily air quality standards of 50 µg m −3 for PM10 was exceeded. The French alert daily threshold of 80 µg m −3 has also been exceeded in the French Alps. At the Passy station in the Arve Valley, PM10 concentrations reached 133 and 104 µg m −3 respectively on 13, 14 December.

Observations
Several set of observations are available to evaluate the model outputs for the different numerical experiments. The 2 m temperature (TMP), relative humidity (RH) and the wind speed (SPD) are used on an hourly basis at the stations in Table 3. Additional meteorological data issued from the Arve Valley provide insightful information on temperature vertical profiles [51,52] to identify thermal inversions. In short, 9, 10 and 11 temperature sensors have been installed along the slope close to Marnaz, Passy and Chamonix, respectively. For the Marnaz and Passy transects, sensors are installed at 3 m above ground level while for Chamonix they are installed on ski lifts pylons (at about 25 m a.g.l.). For the air pollutants a full list of stations from the French national network is available to assess the model on criteria pollutants (PM10, PM2.5, NO, NO 2 , O 3 , SO 2 ). The chemical composition of the PM10 fraction is available on a daily basis for a limited number of stations (Passy, Chamonix, Marnaz and Grenoble) from the DECOMBIO program [53,54] and the CARA program. The dust concentrations are estimated multiplying by a factor of 8 the Calcium concentrations and the Organic Matter by applying a factor of 1.6 from Organic Carbon concentrations [55]. The list of stations is presented in Table 4. All sites are urban, periurban or rural background stations. Usual errors statistics used for the analysis are computed by formula reported in Appendix A.

Numerical Experiments
The reference simulation REF01 is performed over all domains and two other experiments are performed over some selected domains and shorter periods. The description and naming of simulations are provided in Table 5. For the reference simulation the direct and indirect effect of aerosols are turned off (aer_opt = 0). In order to provide to the air quality community some guidance to improve simulations in mountainous areas, five additional simulations are performed. They are designed to evaluate the on-line coupling between air pollutant concentrations and the meteorology through direct and indirect effects (CPL01, CPL02). New WRF parameterization schemes (WRF01, WRF02 and WRF03) are also compared to those of REF01 in order to estimate a possible gain with a focus on near ground meteorological variables and concentrations. All these scenarios for the studied period are initialized and driven at the boundaries (EUR01) with the REF01 simulation.
For the CPL01 simulation, a full on-line coupling is activated. The choice for the coupling approach was to be the least intrusive as possible in the two models and to implement a dialogue between them by using an external coupler, the Ocean Atmosphere Sea Ice Soil-Model Coupling Toolkit (OASIS-MCT) ( [56]). For the direct effects, CHIMERE sends to WRF the Aerosol Optical Depth, Single-Scattering Albedo (calculated with the Fast-JX on-line model) and the asymmetry factor. For the indirect effects, CHIMERE sends to WRF the aerosol number size distribution, the hygroscopic aerosol number size distribution, the aerosol bulk hygroscopicity and the ice nuclei. The variables are exchanged with a frequency of 20 min. This simulation is complemented by a scenario approach as used in [25,57] to better investigate and isolate the role of PM on radiative effects. This experiment, called CPL02, results in switching-off the PM emission from the residential sector for both domains ALP0033 and GRE0011. This activity sector has been selected as it is targeted of major importance in this region in wintertime but other sectors could have been chosen for this test.
The WRF01 simulation is characterized by the switch off of the horizontal diffusion (diff_opt = 0), the change of the scheme for the PBL from the YSU to the Mellor-Yamada-Janjic TKE scheme, and for the surface layer from the revised MM5 Monin-Obukhov scheme to the Janjic-Monin-Obukhov scheme. In recent studies [15,16], some improvements of the WRF simulations were reported just by changing the diff_opt option commonly set to 2 or 4 in the WRF namelist file. The WRF model handles subgrid-scale mixing in the vertical via the PBL scheme and in the horizontal via the diff_opt namelist option. Diffusion can be computed along model surfaces (diff_opt = 1), the default for real-data cases, or in physical space (diff_opt = 2). In the present situation, computing diffusion in physical space is desirable but can require rather high resolution to avoid instabilities associated with sharp terrain gradients.
The default option, however, causes mixing up and down mountain slopes, which has been shown to be problematic in simulations of topographically confined cold pools. This is because adjacent grid points can vary dramatically in height, which in terms of cold pools means sizable differences in temperature, condensate and water vapor. In their study [15,16], they find that deactivating horizontal diffusion entirely (i.e., diff_opt = 0) results in a substantial 11% increase in overnight relative humidity with lower temperature closer to observations. In addition, according [40] and reference therein the YSU PBL scheme is generally characterized by warmer and drier daytime PBLs, which are more consistent with observations, while the MYJ scheme's purely local treatment of larger-scale eddies prevents the PBL from mixing as deeply to produce cooler and moister conditions. Therefore MYJ [58] as a local scheme is chosen in this simulation as it could be more indicated to better capture the cold conditions close to the ground (inversion layers in deep Valleys).
In order to evaluate only the effect of the diff_opt option, WRF02 simulation is a test case only accounting for the change diff_opt = 0. According to the WRF documentation and usual best practices (https://www.climatescience.org.au/sites/default/files/WRF_best-practices.pdf), topographic effects managed by slope_rad and topo_shading options are expected to be important for resolution below 1 km (2 km for slope_rad).
As the target domains have a 1.2 km resolution, a sensitivity test is proposed for the WRF03 simulations for the ALP0033 and ARV0011 domains. The slope_rad option allows for the modification of surface solar radiation flux according to the terrain slope while topo_shading will allow for the shadowing of neighboring grid cells over a maximum distance of 25 km in our case.

Evaluation of the Meteorology
Gridded data from the E-OBS v20.0e dataset [59] are used to evaluate the quality of the meteorological simulation for the 2 m temperature (T2). This dataset is provided at 0.1 • resolution over Europe. Only EUR01 and ALP0033 are evaluated and compared over the ALP0033 area to estimate the added value of a higher resolution. The E-OBS dataset is remapped with NCO [60] over the simulation domains to perform the grid to grid error statistics.
Error statistics displayed in Figure 2 show a general underestimation of the mean temperature with a negative bias of −1.05 K and −0.82 K with spatial correlations of 0.75 and 0.70 on average for EUR01 and ALP0033 domains, respectively. A better bias drives an improved RMSE between APL0033 and EUR01 which are respectively 0.034 and 0.048 K for these domains. The spatial correlation is in general lower for the ALP0033 domain due to the use of a rough NCO bilinear interpolation which does not take into account the altitude. The root mean square error is on average better for the EUR01 domain (not shown here). These negative biases have been reported in previous works [61][62][63]. This result is not surprising since EUR01 and E-OBS have the same resolutions and the EUR01 simulation is nudged with global model fields which account for assimilated observations. In Supplementary Material Figure S2, the chart displays the evolution of the average precipitation simulated over ALP0033 versus E-OBS and clearly the model is able to capture the three rainfall events during the period.
ALP0033 is a free run only constrained by boundary conditions and can locally (particularly over this mountainous area) be different from a smoother fields issued from E-OBS. For both ALP0033 and EUR01 the maximum temperatures are more underestimated than the minimum temperatures. The evolution of the error statistics shows a larger negative bias in November which turns slightly positive during EPII but with a deterioration of spatial correlations. For both domains a negative correlation is observed (about −0.60) between the bias and the correlation, the spatial correlation deteriorates when the bias increases.
Over the Alps, in valleys, the hourly timeseries in Figure 3 confirm the overestimation of minimum temperature, the model is not able to capture the cold pool in Grenoble and Annecy during EPII. In Grenoble the evolution of the 10 m wind speed is rather well reproduced with an improved correlation and RMSE by increasing the spatial resolution. For the relative humidity the error statistics are poor but largely improved with finer resolutions.  Tables 6 and 7 provide an overall evaluation with consistent datasets respectively for the two target domains ARV0011 and GRE0011 for the reference simulation REF01. In order to characterize the effect of the resolution, stations used are those present in each domain only. Results show that a finer resolution does not necessary provide the best error statistics for the PM but for NO 2 there is an obvious added value particularly in the Arve Valley. Nitrate concentrations are on average overestimated particularly at high resolution while sulfate concentrations are underestimated but clearly improved when increasing the resolution. EC and OM benefit from an increase of resolution but OM remains largely underestimated particularly in the Arve Valley. As expected an increased resolution provides higher concentrations close to the sources. A total integrated mass M for a given period T for each species represented by a concentration c(t, x, y, z) in mass per volume at a given time t in a grid cell dx × dy × dz (with dx, dy and dz are respectively the grid dimensions in longitude, latitude and height), is defined as Equation (1):

Chemical Compounds Concentrations
These simulations at different resolutions are an opportunity to test mass budget consistencies. For chemically inert compounds like Elemental Carbon (EC) the total integrated mass for the target domains should provide similar values. As presented in Table 8, over the Grenoble area this budget slightly differs by +4.9% and +3.2% respectively from the 1 ⁄10 • simulation (EUR01) to 1 ⁄30 • (ALP0033) and 1 ⁄30 • (ALP0033) to 1 ⁄90 • (GRE0011). This difference is weak and can be explained (i) by the meteorological fields which are different for each resolution, (ii) the impact the input/output fluxes at the boundaries and (iii) slight interpolation errors for the coarse simulations interpolated over the fine grid. For the nitrate concentration, the difference is also weak, respectively −0.2% and +4.6% from the 1 ⁄10 • simulation (EUR01) to 1 ⁄30 • (ALP0033) and 1 ⁄30 • (ALP0033) to 1 ⁄90 • (GRE0011). These differences are low (and close to the one calculated for a primary species like EC) for such a species involved in a non linear chemistry with certainly a chemical production and loss within the domain not too much impacted by the spatial resolution. The same behavior of the model is observed for the target domain ARV0011 with higher differences between resolutions. Sodium and Dust concentrations are low in this region. The order of magnitude is correctly reproduced for dust with an improved correlation at a higher resolution while the sodium concentrations are underestimated. Ozone concentrations is not a wintertime issue but this species plays an important role on the atmospheric chemistry. As shown in Figure 4 the background concentrations in the altitude site "Le Casset" are rather well captured with a very weak difference between resolutions. For urban sites located in the valleys, the model strongly overestimates the concentrations whatever the resolution, and the correlations are improved when increasing the resolution from EUR01 to ALP0033.   Figure S4 shows the overestimation of nitrate concentrations during EPII and elemental carbon for Passy and Grenoble. Both over the Grenoble area and the Arve Valley, a clear decrease of simulated primary concentrations at the end of EPII (17)(18)(19) is predicted by CHIMERE which is not confirmed by the observations. When looking at criteria pollutants in Figure  S5 the agreement is better with the same discrepancy at the end of EPII.

Timeseries in
As shown in Figure S6, OA is mainly composed of POA particularly over the Alps. For the highest OA concentrations the POA/OM ratio can reach 90%. OM from wood burning (secondary-aged and primary sources) represents a large fraction of OM, over the European domain it can contribute in winter to 50-70% and between 70-90% over the Alps domains Figure S7. These ratios are in line with those reported in the literature for Alpine regions [64] using 14 C methods, while other works like [65] report lower ratio of 68% with a Chemical Mass Balance method applied in Grenoble. Past works showed a contribution of 47-68% at rural areas but excluding the secondary fraction [66].

Impact of the On-Line Coupling
In the reference simulation REF01, no direct effects of aerosols on the radiation calculations are activated. In the CPL01 sensitivity test, WRF uses the aerosol optical properties and particle number distributions from CHIMERE at each physical time step so that the meteorological variables will be affected through a modification of the radiative budget. For the indirect effects, the Thompson and Eidhammer [67] cloud microphysics scheme was replaced by the Abdul-Razzak and Ghan [68] parameterization because CHIMERE uses a sectional treatment. The CPL01 case is run over a short period EPI end of November. EPI is a short PM pollution episode ending by light precipitations (mixed rain and snow) observed in Grenoble (1-2 mm day −1 ) associated to low level clouds largely influenced by an intrusion of humidity from the South-East. Over Grenoble, the impact of the on-line coupling is generally low but a finer resolution increases the effect ( Figure 6). Table 9 reports the average temperature (TMP), the cumulated cloud liquid water content (LWC) and the average vertical diffusion coefficient (KZZ) over the domain GRE0011 for the reference simulation REF01 and the two sensitivity tests CPL01 and CPL02. During EPI, on 29 November, some fog and mist are observed at night in the Grenoble area. For both simulations the resolution has a strong impact on the mean LWC which is about 1.5 times lower at the finer resolution of 1 ⁄90 • compared to 1 ⁄30 • where PM10 concentrations exceed 10 µg m −3 . The mean temperature is less affected. Regarding the effect of the on-line coupling, there is an increase of 5-10% of the LWC when the on-line coupling is activated (CPL01 vs. REF01) mainly due to a cooler temperature. This is the so-called Twomey effect for warm clouds [69,70]. Table 9. Impact of the resolution and the on-line coupling on the average temperature (TMP), the cumulated liquid water content over the domain (LWC) (as computed by Equation (1)) and the vertical diffusion coefficient (KZZ) in 3D on 29/11 for the target domain GRE0011. Two conditions for the averaging process based on the threshold of 10 µg m −3 are applied on PM10 based on the REF01 simulation.  Table 9. Cont. An increase of CCN will not affect too much the LWC but the size and the number of cloud droplets enhancing the "cloud albedo" effect [67]. The influence of the on-line coupling on the radiative budget leads to slightly lower vertical diffusion (10% lower) where the PM10 concentrations are the highest (threshold of 10 µg m −3 ) within the domain, i.e., in the valley bottom while elsewhere the temperatures are slightly higher with an unchanged diffusion coefficient. This phenomenon can be observed in Figure 7. The modification of the radiation budget due to the activation of the aerosol direct and indirect effects with a modification of cloud structures is observed at the bottom of the valley lowering the temperature with elsewhere in the domain and particularly in altitude sites sometimes higher temperatures. This rather important effect over remote sites far from emission sources and polluted sources has been mentioned in recent works [25,72,73] in West Africa. In their works, it results on a subtle balance of physical processes explained by both direct and indirect effects in relation to the presence of the Atlantic Ocean. For our case it is rather different, the comparison of CPL02 with CPL01 helps to understand the processes at work. When emissions in the bottom of the valleys are switched-off the sequence of phenomena are summarized hereafter:

Variable Simulation Condition
Due to a change of concentrations and the direct effect on the radiative budget, temperature and wind speed are slightly modified at the bottom of the valley but sufficiently to initiate a perturbation of wind regimes in the vicinity. (ii) The change of the integrated PM2.5 concentrations over the domain between CPL02 and CPL01 ranges on average from 5.5 % to 0.3 % for "Plain" and "Altitude" sites, respectively. Altitude refers to surface pressures below 750 hPa and Plain for surface pressures above 900 hPa.
The gradient of water vapor mixing ratio (QVAPOR variable in WRF) is rather high in steep areas, then the total column of water vapor can strongly fluctuate from cell to cell in the model. (iv) Water vapor as a radiative forcer contributes significantly to the greenhouse effect, between 35% and 65% for clear sky conditions and between 65% and 85% for a cloudy day. Water vapor concentration fluctuates regionally and locally as shown in Figure S8. [74]. At night the long-wave radiation is one of the most important variable governing the radiative budget, a change of water mixing ratio initiated in the bottom of the valley by small motions immediately modifies the radiative balance.
Over altitude areas nearby the valley, the albedo is generally high about 0.6 compared to 0.25 in the valley. A higher albedo amplifies this effect on the radiative budget with a higher long-wave radiation emitted from the ground, involving modifications of temperature and then enhancing the modifications of local wind regimes and temperatures.
(vi) Under anticyclonic conditions, when the sun rises, the short-wave radiation becomes dominant in the radiative budget and the perturbations due to long-wave radiation slightly disappears but can reappear after the sunset.
Therefore, according to our findings, a change of PM concentrations in the valley can have more impact at night particularly nearby altitude areas (Surface Pressure < 750 hPa) compared to the plain (Surface Pressure > 900 hPa). The change of PM concentrations is the initial cause but the driving process maintaining the perturbation is due to the spatial variability of the water vapor involving a change of temperatures and wind regimes over the steep areas. This local variability is highlighted in Figure S8 showing the MODIS precipitable water vapor layer as the total amount of water vapor in a 5 km by 5 km column of the atmosphere. Absolute differences for some variables are calculated as |CPL01 − CPL02|/|CPL02| expessed in % (except the temperature |CPL01 − CPL02| in • C), the full explanation is provided in Appendix B.
As shown in Figure 8 there is an obvious direct correlation between the perturbations of ground temperature, the upwelling longwave radiation at the top of the model and horizontal wind speed with a coefficient of determination (R 2 ) exceeding 0.8 while the perturbation on PM2.5 is poorly correlated (R 2 < 0.3) with these variables. This support our analysis that PM are only the trigger but not the factor maintaining the perturbations of meteorological variables. There is no specific sign (positive or negative) for the perturbation of wind speed, temperature and water vapor mixing ratio over altitude areas. For the 28-29 November period a diurnal cycle is observed in the model with lower perturbations in the afternoon compared to the night. This diurnal cycle can be visualized in 3D in Supplementary Material Figure S9 for the absolute perturbations of water mixing ratio and vertical wind speed. This kind of "indirect" effect of PM on water vapor mixing ratio is close to what was reported in a recent study [75] at a larger scale showing the effect of Black Carbon on the Asian monsoon over the Tibetan Plateau. Other recent contributions highlight the potential role of water vapor mixing ratio on temperature [76]. Figure 9 shows that for the vertical wind speed, the perturbation occurs over the most complex orographic area of the domain GRE0011 and can extend rather high in altitude, the difference can be positive or negative with dipoles appearing as a consequence of mass and momentum conservation. As a consequence, the PBL height can vary up to 10% on average at night in the plains and up to more than 20% to 30% over altitude areas. In our analysis, the liquid water content and the presence of slight rains are neglected because of the meteorological conditions offering clear sky conditions most of the time. As previously mentioned on 30 November some slight precipitations are observed due to a weak low pressure system developing over the Mediterranean basin. The total rain mixing ratio integrated in time and space in the GRE0011 domain is multiplied by a factor of 8 between the the reference case REF01 and the CLP01 simulation (respectively from 110 to 850 ktons.day). While no precipitation is simulated in the reference simulation, a signal of precipitation even low is visible at Grenoble when activating the on-line coupling in CPL01, Figure 10. The 3D shade shows more precipitations occurring in the morning on the Western part of the domain at the lowest levels where the PM concentrations are the highest the day before. An increase of precipitations is not straightforward to explain because several processes can act in different ways. First, (i) it has been reported in previous studies that producing smaller cloud droplets can delay the precipitations [77], second, (ii) an increase of LWC will mechanically increase the precipitation rates and (iii), an increase of precipitations will scavenge more the particles (at the origin of CCN). Therefore, during a wintertime pollution episode in Alpine regions the aerosol-cloud-radiation interactions can change the quality and quantity of low level clouds and then activate more precipitations in the model with non-negligible impacts on meteorology and air pollutant concentrations through a modification of the dispersion and the deposition processes. CPL01 is supposed to be more realistic and the CPL02 simulation allows to assess the strength of this interaction through the contribution of the residential sector driven by wood burning on meteorological variables. CPL02 provides on average close results to the CPL01 simulation but there are visible differences mainly occurring in the morning over altitude areas with difference on temperature reaching absolute differences close to ±1 • C.

Impact of WRF Configurations
As shown in the section presenting the evaluation of the REF01 simulation, there is a clear overestimation of nighttime temperatures during EPII while for EPI the model behaves better. The vertical temperature profiles in Passy are presented in Figure 11 for two dates and confirm the model discrepancies at the ground level. To compare the model with the observations along the slope, a simple bilinear interpolation of the 2 m temperature is performed. The strong observed thermal inversion in the morning is not reproduced by WRF, a δT = −5 • C is observed between the Passy station and 200 m above while the model simulates a rather flat profile. The WRF02 simulation (only turn diff_opt option to 0) seems to provide the lowest bias. These results for two dates are reprsentative of the dicrepancies of the model which is not able to represent the strong thermal inversion. For the highest resolution, the effect of different WRF settings is compared to available observations for the Annecy station (LFLP). As shown in Figure 12, there is no clear advantage of one setting compared to others. The WRF03 simulation setting was supposed to better account for topography effect but does not provide an improvement of time correlation on hourly temperatures which are largely overestimated at night for all configurations, however the biases are improved for the WRF03 simulation for all variables. No PBL observations are available in Annecy but we can see large differences between WRF configurations. The MYJ PBL scheme produces on average the lowest PBL height (PBLH) (about 59 m) with a minimum PBL higher than the YSU scheme (visible during the nights). During this episode the wind speed is overestimated by a factor of 2 with a poor correlation whatever the configuration. However, for such low values, the rounding process of observations (potentially associated to unit change) could amplify this effect, the observation data processing could at least explain on average 0.5 m s −1 or knot (depending on the original unit) of the positive bias. Switching off the diff_opt option (WRF02) only compared to REF01 has a rather important impact on wind speed but is not necessarily better. diff_opt=0 in combination with diff_6th_opt=0 does not allow horizontal diffusion along the slopes and this option should produce substantially more fog development at the ground level [16] but unfortunately in the studied period no significant fog is recorded. The WRF02 case provides the highest wind speed and then the highest boundary layer heights.
Average values of some meteorological variables for the various configurations are reported in Table 10 for the most polluted day during episode EPII. Lower temperatures seem associated to higher LWC.
The change of meteorological simulations largely affects the computation of concentrations. WRF01 and WRF02 improve the bias for nitrate and EC concentrations which were overestimated in the REF01 simulation for the Passy station located in the Arve Valley Figure 13. Surprisingly, the WRF01 produces the lowest concentrations with the lowest KZZ and PBL values. Actually, as previously mentioned the YSU scheme is capped by a minimum PBL value higher at night allowing lower nocturnal concentrations which is confirmed by the hourly timeseries of NO 2 and PM10. For all configurations, sulfate remains underestimated, this species in the model is mainly influenced by long range transport and certainly local sources are not taken into account due to gaps in the emission inventory. The WRF03 simulation is supposed to better account for topographic effects and the main effect results in an increase of concentrations which is rather fair for OM and NO 2 but deteriorates the bias for Nitrate and EC. For the Marnaz and Chamonix stations results are similar with a less pronounced impact on nitrate concentrations (Supplementary Material Figures S10 and S11). Figure 13. Impact of WRF settings on species concentrations at 1 ⁄90 • resolution (ARV0011) for the EPII period at Passy. Unfortunately the large positive bias on nitrate is not fully improved by one of the WRF configuration. As NO 2 concentrations are fairly reproduced a too large NOx emission is not the main reason of this discrepancy on nitrate concentrations. Looking at the daily time series of nitrate concentrations for Grenoble, Passy, Marnaz and Chamonix for the full period the large positive bias during EPII in December which is not observed in November seems correlated with the overestimation of temperatures. Therefore, the chemical production of nitric acid could be the main reason due to its sensitivity to meteorological conditions.

Transport of Air Pollutants
Evaluating and analyzing the transport of pollutants in an environment with such a complex topography is not straightforward. This work proposes to analyze the transport of EC, OM and Dust air pollutants through the deposition processes. These species are of major importance for health concerns and to understand the evolution of glaciers melting in mountainous regions. Dust is emitted from remote arid areas, the Sahara contributing in Europe during outbreaks or as a background disseminated in the atmosphere over the North Hemisphere. In the previous sections we showed that the order of magnitude of Dust concentrations was rather well captured by the model. EC and OM are influenced by local sources (particularly EC) in winter and could affect nearby altitude sites, however the simulated PBL heights are rather low during the pollution episodes and the observed low thermal inversion layers should keep the pollutants in the lowest layers of the atmosphere. Figure 14 shows the accumulated deposition over the ARV0011 domain from 15 November to 20 December. Clearly, the wet deposition of desert dust dominates over the glaciated area of the Mont-Blanc massif, the amount of 40 mg m −2 is coherent with what was usually observed in this area [78]. The deposition of EC does not exceed 1 mg m −2 and is in line with what is commonly observed over other mountainous areas [79,80]. EC and OM deposition are the highest at the West and East edges of the Glacier. The 6-hourly timeseries of wet and dry deposition for EC, OM and Dust over the Glaciers surrounding the Mont-Blanc is provided in Supplementary Material Figure S12 and can explain the dynamic of the air pollutants advection.
For the two main synoptic rainy periods around the 20 November and 19 December, EC dry deposition occurs before or during the precipitation event. EC concentrations is less affected by wet deposition during large scale rainy periods because the sub-cloud scavenging is less effective than for coarse Dust particles since EC is prevalent in fine fraction of particles. In addition, under these meteorological conditions EC emissions are not collocated with fog or clouds. The synoptic flux blows-up the accumulated EC in the valley during the PM pollution episodes and transport the pollutants up to the Glaciers and the first areas affected are the West or East edges of the Glacier. The 30 November event is different, the precipitations are lower with a lower synoptic flux and the deposition over the Glacier is mainly driven by wet deposition of EC advected from the surroundings valleys within cloud droplets.

Conclusions
This work aims at evaluating the capacity of a chemistry transport model like CHIMERE to simulate the air quality in a complex orographic area. The French Alps are often submitted to pollution episodes particularly in wintertime with a large contribution of wood burning. Our findings can be summarized hereafter.
CHIMERE is able to capture the evolution of the PM and NO 2 concentrations even if some constitutive species like nitrate are overestimated at local scale while the organic matter remains underestimated. In some cases WRF is not able to retrieve the thermal inversions observed in December 2013, despite the use of other WRF settings expected to improve the results as reported in the literature. The use of different WRF configurations greatly affects the air pollutant concentrations with differences which can reach up to 30-50% for some PM constitutive species. Finally the reference simulation REF01 with the original CHIMERE and WRF configurations provides satisfactory results.
For the Arve Valley, an increase of the resolution is certainly desirable but difficult to support for operational uses. There is an added value of an improvement of the horizontal resolution. Increasing the resolution affects the spatial distribution of pollutants increasing the concentrations of primary pollutants over emission zones and decreasing the concentrations in remote areas. In terms of integrated mass there is a rather good conservation even for secondary pollutants where a maximum of 10% differences are observed between the finer and the coarser resolution.
The on-line coupling is expected to be closer to the reality. It slightly affects the simulation of concentrations and the meteorology. Statistically there is no large difference when activating the on-line coupling. However, thanks to the coupling, precipitations are enhanced in our case and the model results match better with the observations, even if the intensity remains too low in the model. Surprisingly our finding suggests a non-negligible impact of the presence of PM in the valley with the meteorology in the surrounding altitude areas. PM concentrations are the trigger of this phenomenon through the direct and radiative effects, but the motion of water vapor is the main driver of these perturbations that operates at night, and disappear during daytime. The sharp gradients of water vapor in mountainous regions drive and amplify the perturbations over high altitude sites.
In wintertime, the altitude sites can be affected by the pollution in the valley. PM and particularly EC can be deposited over the Glaciers through the dry deposition fluxes first and by precipitations.
Over the Alps, the wet deposition of background desert dust is of major importance compared to EC. However the radiation absorption efficiency of EC is more important with an imaginary part of the complex refractive index of 0.569 compared to 0.008 for Dust.
This study showing how an initial small perturbation (PM emissions) through the radiative effects can be propagated and amplified in a mesoscale meteorological model is questioning and raises the issue of the "hyper" sensitivity of a model to a parameterization and possible associated divergences in case of free run not constrained by nudging processes and/or boundary conditions. In addition, as a perspective, an improvement of the resolution of the CHIMERE-WRF suite could help to solve the problem of thermal inversion layers that are sometimes poorly captured by the model. Data assimilation of local meteorology could help in such deep valleys to better constrain the model.

Supplementary Materials:
The following are available on-line at http://www.mdpi.com/2073-4433/11/6/565/ s1, Figure S1: Synotic meteorological situation November-December 2013 for some key dates from metoffice.gov.uk, Figure S2: Evolution of the average precipitation simulated for the REF01 simulation over ALP0033 comaped to E-OBS observations, Figure S3: Impact of the horizontal resolution 1/10 o -1/30 o -1/90 o (from left to right) on daily concentration fields on 11/12 for six pollutants and species (PM10, NO 2 , OM, EC, Nitrate, Sulfate from top to bottom) over the Grenoble area (GRE0011). Colored squares represent the observations, Figure S4: Daily time series of chemical species for the full period in Grenoble, Passy and Chamonix, Figure S5: Daily time series of chemical species for the full period for PM10, PM2.5 and NO 2 for key stations over the ALP0033 domain, Figure  S6: Mean simulated POA/OM ratio by minimum threshold on OM averaged over each domain. The encapsulated chart is the evolution of the average ratio over the studied period, Figure S7: Evolution of the fraction of wood buring OM in the total OM fraction over the studied period, Figure S8: MODIS Precipitable Water Vapor layer shows the total amount of water vapor in a 5 km by 5 km column of the atmosphere, measured in centimeters (cm). Color scale ranges from blue to yellow (0 to 1.4 cm), Figure S9: 6-Hourly evolution on 29/11/2013 of the difference CPL01 -CPL02 for the water vapor mixing ratio in kg kg −1 (on the left) and the vertical wind speed in m s −1 (on the right side), Figure S10: Impact of WRF settings on species concentrations at 1 ⁄90 o resolution (ARV0011) for the EPII period at Marnaz, Figure S11: Impact of WRF settings on species concentrations at 1 ⁄90 o resolution (ARV0011) for the EPII period at Chamonix, Figure S12: 6-hourly time evolution of the wet and and dry deposition for DUST, EC, OM and SOA species over the Glaciers surrounding the Mt-Blanc, Table S1: List and location of the main air quality monitoring stations (PM10, PM2.5, NO, NO 2 , O 3 , SO 2 ) for error statistics, Table S2: List and location of the main air quality monitoring stations (PM10, PM2.5, NO, NO 2 , O 3 , SO 2 ) for error statistics (next), Table S3 Acknowledgments: The authors acknowledge Laurent Létinois (INERIS/LCSQA) for extracting the official air quality data from the GEOD'AIR database (available on request at https://www.lcsqa.org/fr/les-donneesnationales-de-qualite-de-lair). Synoptic meteorological data (MESONET network) is available at http://www. weather.uwyo.edu/surface/meteorogram/europe.shtml. Meteorological data at Grenoble Le Versoud are available at http://www.meteociel.com. We acknowledge the E-OBS dataset from the EU-FP6 project UERRA (http://www.uerra.eu) and the Copernicus Climate Change Service, and the data providers in the ECA&D project (https://www.ecad.eu). The PM10 chemistry and the vertical temperature profiles in the Arve Valley were acquired during the DECOMBIO project funded by Ademe and Primequal (1362C0028) with the dedicated efforts of F Chevier, J Allard, and C Coulaud. The PM10 chemistry data for Grenoble is funded by the Program CARA and internal fundings at IGE. Many thanks to numerous people at IGE for the chemical analyses of the Arve and Grenoble series, and to many people at Atmo AuRA for the collection of these samples. We acknowledge the use of imagery from the NASA Worldview application (https://worldview.earthdata.nasa.gov), part of the NASA Earth Observing System Data and Information System (EOSDIS). For the NASA Earth Observations of MODIS water vapor, special thanks and acknowledgments go to David Herring, Jesse Allen, Jacques Descloitres, Gene Feldman, Chelle Gentemann, Mark Gray, Goran Halusa, Jackie Kendall, Norman Kuring, Rebecca Lindsey, Chris Lynnes, Alex McClung, Bill Ridgway, Holli Riebeek, Jeff Schmaltz, Michon Scott, Rob Simmon, Allan Smith, Reto Stockli, Ross Swick, Stephanie Schollaert Uz, Tom Whittaker and Wenli Yang for generously sharing their time, expertise, imagery, and data. Very special appreciation goes to Vince Salomonson and Michael King for their initial support and encouragement and Steve Platnick, EOS Project Scientist, for his continuing support.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Error Statistics
The mean bias (MB), Pearson correlation (R) and the root mean square error (RMSE) are defined herebelow for n the number of data, M the model predicted value and O the corresponding observation.

Appendix B. Averaging and Cumulating Computation
Cumulated 3D variable V(t) of variable v(t, x, y, z) over x,y,z for a given time t: x,y,z v(t, x, y, z) dx dy dz (A5) Averaged 2D variable V(t) of variable v(t, x, y) over x,y with n the number of grid cells (x, y) for a given time t: A mask on the surface pressure P sur f is applied for the averaging and cumulating process to target Plain and Altitude areas according to: • P sur f < 750 hPa for Altitude areas • P sur f > 900 hPa for Plain areas For most variables V the average differences or perturbations (p in %) between two scenarios indexed (1) and (2) for a given time t are computed as: except for the Temperature in • C (as a relative value is not relevant here) which is computed as: