Atmospheric Pollutant Dispersion over Complex Terrain: Challenges and Needs for Improving Air Quality Measurements and Modeling

Abstract

Pollutant dispersion processes over complex terrain are much more complicated than over flat areas, as they are affected by atmospheric interactions with the orography at different spatial scales. This paper reviews recent findings and progress in this field, focusing on both experimental and modeling perspectives. It highlights open questions and challenges to our capability for better understanding and representing atmospheric processes controlling the fate of pollutants over mountainous areas. In particular, attention is focused on new measurement techniques for the retrieval of spatially distributed turbulence information and air quality parameters, and on challenges for meteorological and dispersion models to reproduce fine-scale processes influenced by the orography. Finally, specific needs in this field are discussed, along with possible directions for future research efforts.

1. Introduction

Many mountains in the world are undergoing a process of increasing anthropization. Indeed, nowadays roughly 900 million people, i.e., 13% of the world’s population, live in mountainous regions [1,2]. An even higher number of persons are attracted to mountainous areas for tourism and recreational activities every year. For example, in the Alps, an overall demand of more than 464 million overnight stays per year is estimated [3]. Moreover, populations living in mountain regions increasingly concentrate in rapidly expanding urban areas; it is estimated that roughly 30% of the world’s mountain population lives nowadays in cities, and this percentage reaches 50% in developed countries [2].

Mountain valleys and passes are also often crossed by main transport routes connecting different countries. For example, the Alps are crossed by 4200 km of main roads, with 6 million vehicles and 1.24 × 10¹¹ kg of freights being transported across the Alpine range every year [4]. As a consequence, the anthropogenic impact on mountainous areas raises concerns from many viewpoints. In particular, one of the most serious hazards involves pollutant emissions and their effects on air quality. Major sources of pollutants in mountainous regions are generally located inside valleys and basins, where the main urban areas and traffic routes develop, and most of the population lives. This is shown, as an example, in Figure 1, displaying annual emissions of NO_x and PM₁₀ in South Tyrol, in the Italian Alps [5]. Under unfavorable meteorological situations, pollutants may accumulate at the bottom of valleys and basins, implying that a large part of the population is potentially exposed to high concentrations. In fact, in such conditions, atmospheric dispersion, i.e., the combined effect of transport/advection by the mean wind and diffusion due to turbulent mixing, may be strongly inhibited. High levels of pollutants are harmful, not only to humans, wildlife and livestock, but may also cause damage to the environment, especially vegetation and forest heritage [6]. Furthermore, pollutants may be transported by orographic flows towards mountainous areas from adjacent regions, where they are emitted. As a consequence, significant concentrations of harmful species, especially secondary pollutants, may also be found in areas where they were not locally emitted [7,8,9].

The relevance of the problem of pollutant dispersion has been a primary driver for many research efforts focused on atmospheric processes over complex terrain in recent decades. In fact, the need for better understanding and suitable modeling of the fate of pollutants in mountainous regions stimulated a number of projects, which involved both experimental activities, such as extensive field campaigns in complex environments, and the advancement of new modeling capabilities over complex terrain, encompassing both meteorological and air quality models. Among the projects primarily focusing on the study of air pollution dispersion over complex terrain, it is worth mentioning Monitoring and minimization of traffic-induced noise and air pollution along major Alpine transport routes (ALPNAP) [10,11,12], Bolzano Tracer Experiment (BTEX) [13,14], Persistent Cold-Air Pool (PCAPS) [15], Transport of Air Pollutants over Complex Terrain (TRACT) [16], Mesoscale transport of atmospheric pollutants across the central Alps (TRANSALP) [17,18], Vertical Ozone Transports in the Alps (VOTALP) [19], and Vertical Transport and Mixing (VTMX) [20]. Other projects aiming at characterizing dispersion in complex terrain and directly or indirectly focusing on air pollution will be reviewed in this paper. Results from these projects have greatly contributed to the advancement of knowledge in this field. However, several aspects, concerning both our understanding and our ability of observing and modeling processes related to pollutant dispersion over complex terrain, still remain to be investigated. In fact, the understanding and the quantification of pollutant dispersion over complex terrain are much more difficult than over flat areas, as dispersion processes are affected by atmospheric interactions with the orography at different spatial scales. In particular, atmospheric flows over complex terrain are characterized by a continuous and interacting range of scales, from synoptic forcing to mesoscale circulations and turbulence fluctuations. In complex terrain, the mechanical and thermal influence of the orography can modify the large-scale flow and produce smaller-scale motions which would not exist on flat terrain [21,22,23], thus enhancing the spatial and temporal variability of atmospheric processes relevant for pollutant dispersion.

Based on these considerations, this paper aims at reviewing recent research activities focusing on several aspects connected to pollutant dispersion over complex terrain, from both experimental and numerical modeling perspectives. A review of recent progress in the field will lead to open questions and challenges for improving our capability to better understand and represent atmospheric processes responsible for the fate of pollutants over complex terrain. Following this objective, the present paper is organized as follows. Section 2 contains a general overview of the processes affecting pollutant transport and diffusion over complex terrain, focusing in particular on the peculiarities of atmospheric processes over mountains when compared to flat terrain. Section 3 is dedicated to measurement techniques supporting air quality management and the understanding and modeling of atmospheric dispersion processes, with a focus on new observational techniques, remote sensing, and innovative low-cost sensors for air quality monitoring. Section 4 and Section 5 concentrate on numerical modeling and provide an overview of the state of the art and perspectives of both meteorological and dispersion modeling. Finally, conclusions and an outlook for future research needs are given in Section 6.

2. Atmospheric Processes and Pollutant Dispersion over Complex Terrain

The concentration of pollutants in the atmospheric boundary layer (ABL) depends not only on emissions, but also, and sometimes crucially, on atmospheric dispersion. High concentrations of pollutants therefore follow not only from strong emissions, but also from reduced dispersion, especially in the atmospheric layer adjacent to the surface, where most of the pollutants are emitted.

2.1. Temperature Stratification

Temperature stratification (static stability) is one of the main factors controlling turbulent mixing in the atmosphere: the most critical situations occur under wintertime anticyclonic episodes, when radiative cooling at the ground leads to the formation of a very stable air layer above the surface. This situation is usually referred to as ground-based inversion; it inhibits turbulent motions, and consequently keeps pollutants “trapped” in that layer [24,25,26,27,28,29]. This situation may last for several days to several weeks, and is more frequent in mountainous regions; here the reduced sky-view factor limits the penetration of solar radiation into valleys and basins, and the nocturnal drainage winds favor the accumulation of cold air at low levels.

2.2. Mountain Boundary Layer Height

A key parameter affecting pollutant concentration is the height of the ABL, as it determines the volume of atmosphere available for pollutant dispersion. Turbulent mixing comes again into play, since it also contributes to controlling the evolution of the ABL depth. Over flat terrain, the depth of the ABL can often be assumed to be horizontally homogeneous, whereas over complex terrain the so-called “mountain boundary layer” (MoBL) [30] is spatially highly inhomogeneous and interacts with mesoscale flows. As a consequence, an accurate determination of the height of the MoBL is a more difficult task than over flat areas [30]. In fact, apart from quite substantial uncertainty arising from data quality and parameter choices even over flat terrain (cf., e.g., [31]), at least well-established criteria are available in the literature. Conversely, more interacting factors contribute to the spatio-temporal variability of the MoBL height in complex terrain, as extensively discussed in a recent review paper [32]. For example, heterogeneity in the penetration of solar radiation in mountain valleys and basins can produce significant spatial variability in the surface energy budget, influencing both atmospheric stratification and the boundary layer height [33,34]. Additionally, quantities defining the reference “boundary layer height” based on classical (i.e., flat terrain) concepts under convective conditions, do not, in general, correspond to the height up to which mixing is effective (i.e., the MoBL height). For example, Lang et al. [35] highlighted, using idealized large eddy simulations (LES), that different definitions of the boundary layer height yield similar results over flat terrain (Figure 2a), while significant discrepancies arise when idealized topographic features are introduced into the simulation (Figure 2b). Likewise, under stable stratification, the height of the boundary (inversion) layer can be modulated by the interaction of the upper-level wind with down-valley flows of tributary valleys, as the direction of the upper-level wind varies during the stable episode [36].

2.3. Strong Synoptic Forcing

When the upper-level wind speed is considered, two main situations affecting pollutant dispersion in the lowest layers in mountainous regions can be distinguished. One situation occurs with strong upper-level wind speed, associated with a significant synoptic forcing. Over a plain, these conditions usually favor the dispersion of pollutants. However, in mountainous terrain, depending on the alignment of the upper-level wind with the valley axis [38], the valley bottom may be sheltered from downward intrusions of the wind from above. Additionally, the penetration of the upper-level wind down to the lowest levels of valleys and basins may be further prevented by underlying pre-existing and persistent stable layers, such as those associated with nighttime cooling, especially in wintertime [39]. Under these conditions, the flow inside valleys is decoupled from the airflow situation above the crest level, and the persistence of such conditions (i.e., light local flows and stable stratification) favors pollutant accumulation at the lowest levels.

2.4. Weak Synoptic Forcing and Thermally Driven Flows

The second situation occurs under persistent anticyclonic pressure fields, generally associated with weak synoptic forcing and clear skies. In this case, the circulation inside the valley or basin is driven by the local heating and cooling of the air layers adjacent to the ground, leading to down-slope and down-valley flows during the night and up-slope and up-valley flows during the day [40,41]. The slope and valley flows control the transport and the diffusion of pollutants. In summer, the up-slope and up-valley flows are generally well developed, and may be associated with convective cells, which can extend very high above the valley bottom, up to a few kilometers, thus transporting pollutants to high altitudes [42]. In this case, up-slope circulations strongly enhance the mass exchange between the valley boundary layer and the free atmosphere with respect to pure turbulent exchange [43]. On the contrary, in wintertime, due to the reduced heating of the valley atmosphere, down-slope and down-valley flows dominate the valley circulation [44,45]; as shown in [46] from numerical simulations in the Grenoble Valley, this “nighttime” flow regime may last about 18 h (over 24 h), ceasing only for a few hours around midday, during which a weak convective activity may be observed above the valley bottom over a height of, at most, 50 m. In these situations, wind speed is generally weak (less than 1.5–2 m s⁻¹), with frequent calm regimes (with wind speed less than 0.5 ms⁻¹) and an inversion layer characterized by a positive temperature gradient can form in the valley up to an altitude close to the crest level. As a result, pollutant dispersion is strongly reduced during a wintertime episode [47,48]. In low wind speed conditions, pollutant dispersion is governed by meandering, a phenomenon still largely under study (see a brief review in [49]), dispersing plumes over wide angular sectors [23].

The inversion layer, typical of these wintertime episodes, along with topographic effects, leads to the formation of stagnation and ventilation zones [50,51]. In the sections of the valley displaying a nearly basin-like structure, stagnation zones develop in which pollutants travel over short distances from their emission points [11,52]. In these stagnation zones, pollutant concentrations increase steadily in time during the episode [53,54]. On the other hand, in ventilation zones, such as a valley section opening onto a plain, ventilation effects lead to the evacuation of pollutants [29]. The end of the persistent anticyclonic wintertime episode occurs when the upper-level wind is able to reach the valley floor and/or the energy input is sufficient to break up the inversion layer [55]. As a result, pollutants can be transported above crest level by vertical fluxes or out of the valley by the channeled wind. An example of the progressive accumulation of pollutants on the valley floor during an anticyclonic wintertime episode is reported in Figure 3. This figure displays a time–height plot of wind speed measured by a SODAR (SOnic Detection And Ranging), along with PM₁₀ concentration, during a field campaign of the ALPNAP project in the Alpine Adige Valley [11]. It can be noted that PM₁₀ concentration remains low in the first three days, when moderate winds reach the valley floor. On the other hand, a progressive accumulation of PM₁₀ occurs in the subsequent days, characterized by weak winds, until a progressive increase in the wind speed terminates the air pollution episode.

Thermally driven circulations affect the transport of pollutants in many other different ways. For example, valley winds may advect highly polluted air downstream from major pollutant sources. Conversely, valley circulations may also have a cleansing effect, when cleaner air is transported into polluted regions. More generally, pollutant transport by a down-valley flow depends both upon the stratification in the valley and upon local orographic features. Indeed, if a strong ground-based inversion exists at the valley bottom, the down-valley wind may detrain at its level of neutral buoyancy, leaving the valley bottom unperturbed. The detrainment may be also associated with the excitation of internal gravity waves [36,56,57], which may result in fluctuations of the pollutant concentration [29]. On the other hand, up-slope or up-valley winds may advect primary pollutants or precursors (also from natural sources, such as forests) to upper levels in the atmosphere; there, the different exposure to stronger radiation and/or other ambient conditions may affect (photo)chemical reactions, producing secondary pollutants (e.g., ozone), which eventually get drained to lower levels by nocturnal down-slope flows. At larger spatial scales, plain-to-mountain circulations can advect pollutants and their precursors towards mountain ranges [58], where they may be embedded into smaller scale circulations, and eventually transported into the free atmosphere (a process referred to as “mountain venting” [59,60,61]), also depending upon the geometry of the mountain range.

The development of thermally driven circulations is strongly affected by land cover heterogeneities, which are more common in mountainous areas. Changes in surface properties (such as albedo, thermal capacity, and conductivity) between valley floors, sidewalls, and mountain tops may modify surface–atmosphere heat exchanges and thermal contrasts driving these local circulations. Other climatic influences being equal, mountainous areas are more likely to have snow-covered terrain for a significant part of the year. Snow cover has many effects on the surface layer. First, snow cover on the ground affects the thermal structure, favoring stable stratification. Moreover, irregular snow cover, typically associated with the different exposure to the solar radiation of the valley sidewalls or with the varying altitude of the valley floor, may also produce asymmetric local circulations, influencing the redistribution of pollutants [52]. The higher amount of reflected radiation from extended snow-covered areas may also favor photochemical reactions, leading to an increased production of secondary pollutants [62,63].

2.5. Urban Areas

Local circulations may also be affected by land cover heterogeneity associated with extended urban areas. Indeed, cities usually experience higher near-surface temperatures, especially during nighttime [64], leading to the so-called “urban heat island” effect. Higher urban temperatures can cause the development of an “urban breeze”, with convergent motions toward the city center at low levels [65]. Thermal hot-spots induced by urban areas may superimpose to the temperature gradients developing inside valleys, and locally modify the normal along-valley horizontal pressure gradients driving the diurnal pattern of up- and down-valley flows [66]. This usually implies a wind convergence over the city, especially during nighttime, when the urban heat island is stronger, thus favoring the accumulation of air pollutants [67]. Moreover, the higher roughness associated with buildings may produce an internal urban boundary layer [68], preventing or limiting the penetration of flows and the ventilation into the urban canopy [69,70]. Similar to the action of upper-level winds on the circulation in a valley, the orientation of urban canyons with respect to the flow direction may determine very different internal circulations, resulting in better ventilation in streets aligned with the mean flow. As a consequence, in cities lying within valleys, streets aligned along the main valley axis are usually better ventilated. Urban areas can also modify the development of the ground-based temperature inversion at night [71,72], as well as its break-up in the morning, with a significant impact both on air quality and on the cross-valley wind system [73,74]. All these aspects still make the case of urban areas over complex terrain a particularly challenging situation.

3. Observational Techniques Supporting Air Quality Management and Pollutant Dispersion Modeling

3.1. Meteorological Measurements

For the assessment (i.e., interpretation, understanding, and modeling) of pollutant dispersion processes anywhere in the world, mean flow characteristics (such as wind speed and direction, stratification, etc.) are required for characterizing advection, as well as information on the turbulence state, which is needed to evaluate diffusion processes. Measurements to assess mean meteorological variables are usually available from meso-networks—albeit generally not at the spatial resolution one would wish. This is particularly true over mountains, where the maintenance and operation of meteorological networks can be highly demanding [75]. On the other hand, measurements of turbulence variables to assess air quality in mountainous terrain are largely unavailable from routine networks. Datasets for these variables therefore stem from research projects with limited duration. One of the first projects was Atmospheric Studies over Complex Terrain (ASCOT) [76], in 1984—even if this study only had a (small) number of “representative” turbulence sites. More recent projects that focused on turbulence and exchange processes in complex terrain—or at least had a strong component on this topic—include ALPNAP in 2005–2008 in the Alps [10,11,12], Boundary Layer Late Afternoon and Sunset Turbulence (BLLAST) in 2011 in southern France [77], Convective and Orographically induced Precipitation Study (COPS) in 2007 in southern Germany [78], Mesoscale Alpine Programme–Riviera (MAP) in 1999 in the Riviera Valley [79], Mountain Terrain Atmospheric Modeling and Observations (MATERHORN) in 2012–2015 in Utah [80], Meteor Crater Experiments (METCRAX and METCRAX II) in 2006–2009 and 2013–2016 in Arizona [81,82], Persistent Cold-Air Pool (PCAPS) in 2010–2011 in Salt Lake Valley [15], Terrain-Induced Rotor Experiment (T-REX) in 2003–2006 in Owens Valley [83], and Vertical Transport and Mixing (VTMX) in 2000 in Salt Lake Valley [20]. All these projects resulted from collaborative efforts and hence were able to measure multiple sites in order to address the problem of spatial variability, which is inherent to boundary layer processes in non-flat terrain [84].

As for long-term measurement installations for turbulence-related variables, there are a handful of “fluxnet towers” (as part of global CO₂ monitoring activities) [85] situated in mountainous terrain (e.g., in Switzerland: Davos [86], in Italy: Ritten/Renon [87], Lavarone [88], and Mt. Bondone [89], and in Austria: Neustift [90]). With the “i-Box” in the Austrian Inn Valley [91] there is, furthermore, a long-term program running, which currently aims at, among others things, providing turbulence data from six sites, which have been selected for their representativeness of particular types of surfaces (valley floor, slope (exposition and slope angle), mountain top, etc.).

Turbulence measurements over complex, mountainous terrain are generally performed with the same instrumentation as over flat terrain, i.e., sonic anemometers/thermometers (point observations), scintillometers (spatially averaged fluxes), and more recently, wind profilers (lidars, sodars) and airborne observations. Due to the terrain complexity, they all exhibit their specific challenges, which are briefly summarized in the following section.

3.1.1. Point Measurements

The post-processing of sonic data includes quality control, coordinate rotation, and the correction of known physical effects (e.g., Webb-Pearman-Leuning (WPL) correction), as summarized and discussed in [22,92,93]. In particular, Falocchi et al. [22] and de Franceschi et al. [93] suggested methods to separate low-frequency components, associated with mean flow and larger vortexes, from smaller scale turbulence responsible for turbulent diffusion. Additionally, over sloping terrain, the coordinate rotation requires particular attention. Both the double rotation (DR) and planar fit (PF) approaches rely on a coordinate system that is normal to the surface [94]. Given the instrument layout of most sonic anemometers (having the least systematic errors for flow parallel to the instrument), these instruments should be installed with the best viewing side normal to the local surface. No clear preference could be found so far concerning DR vs. PF, except that PF seems to have advantages for climatological applications (e.g., long-term CO₂ flux assessment) while DR better corresponds to the assessment of, e.g., local similarity relations. In addition, PF may require more than sectorial planes [92]. For the purpose of comparing measurements to model output, it is relevant that they both refer to the same coordinate framework (vertical vs. normal to the slope), since along-slope turbulent heat fluxes can contribute to the vertical (as opposed to normal) turbulent heat flux considerably, thus influencing stability [95,96].

Many of the “best practice procedures” for the analysis of turbulence measurements in complex terrain are relatively recent and based on the respective authors’ datasets (site characteristics, season, duration, etc.). Therefore, it would be highly desirable to work out an overall best practice, based on systematically comparing different approaches, and their combinations, using different datasets.

3.1.2. Scintillometry

The retrieval algorithm of scintillometry depends on the applicability of scaling relations. Generally, commercial instruments use similarity relations for the temperature structure function based on Monin–Obukhov similarity theory (MOST) (which is strictly valid for the horizontally homogeneous, stationary surface layer only) to derive spatially averaged turbulent fluxes. Thus, the application of scintillometry strongly depends on the applicability of generic MOST relations in complex terrain. Weiss [97] and Weiss et al. [98] found relatively favorable correspondence between scintillometry-derived and point measurements (heat flux) over a gentle local slope. Pianezze [99] reported substantial deviations from “ideal behavior” for the structure functions depending on cross-sections over a steep Alpine valley (the Mesolcina Valley in southern Switzerland). Ward [100] listed a (small) number of additional applications of scintillometry over gentle topography. From Ward’s considerations, it can be concluded that it is particularly important to use an accurate effective beam height, which may be difficult to determine over complex terrain. Hartogensis et al. [101] discussed how to calculate the effective beam height, taking into account varying beam height along the path and its impact on the measured structure parameter (they do not, however, really address any complex orography in this paper).

For a successful application of scintillometry in mountainous terrain, therefore, it is necessary to assess in more detail the applicability (conditions, situations) of MOST relations for the structure parameters, as well as to possibly derive more appropriate parameterizations. Given the potential of scintillometry for model validation (the comparison of inherently spatially averaged turbulent fluxes from model grids with spatially averaged observations), the further development of this area of research is highly desirable.

3.1.3. Lidars

Wind and temperature/humidity lidars are powerful tools to assess the temporal/spatial structure of meteorological variables near the surface, e.g., [102], which is especially useful in complex terrain with its inherently inhomogeneous flows, e.g., [103]. Furthermore, lidar technology allows, in principle (and under suitable conditions, e.g., stationarity during the averaging interval), the retrieval of turbulent (integral) statistics, such as the vertical velocity variance. Depending on the number of instruments involved and the scan strategy, profiles or point information at a specified location in the atmosphere can be retrieved [104]. Retrieval algorithms and the corresponding accuracy [105] of the obtained turbulence information is an area of active research even over flat terrain and has received an increased interest since the availability of commercial wind lidars.

There are only very few applications to date using lidars to retrieve turbulence statistics in complex terrain, e.g., [106,107,108,109]. For example, Adler and Kalthoff [110] have used profiles of vertical velocity moments (variance and skewness) to determine the local MoBL height. Still, the application of lidars to obtain turbulence information in complex terrain is at its infancy and—based on the importance of such information for the simulation or understanding of pollutant dispersion processes in complex terrain—requires a major boost in the near future.

3.1.4. Airborne Observations

Airborne measurements constitute another well-suited platform to detect, not only the mean thermal and dynamical structure, but also turbulence properties controlling the diffusion of pollutants over complex terrain. As to the mean variables, de Franceschi et al. [111] and Laiti et al. [112,113,114] have successfully retrieved the vertical thermal structure associated with valley winds in Alpine valleys by means of measurement flights along carefully designed trajectories, and subsequent data analysis including spatialization via geostatistical techniques (kriging). As to turbulence properties, in principle, traditional methods (Fourier transform of time series, Reynolds decomposition) can be used—with a number of caveats concerning stationarity and homogeneity—to retrieve turbulent statistics. In relation to topographic flows, flights are often conducted relatively far above the terrain, e.g., [115,116] for obvious safety reasons. Weigel et al. [117] have used data from a light research aircraft within and above a relatively small valley to investigate the turbulent kinetic energy (TKE) structure using the Reynolds decomposition technique along flight paths and found favorable agreement with surface-based tower data for “fly-bys”. The wavelet transform has been used to obtain information on the spatial variability of, e.g., turbulent fluxes from airborne data, e.g., [118]. Even if some caveats exist here as well [116], Baur [119] has successfully used this approach to investigate the spatial structure of turbulence characteristics in the Inn Valley located in Austria. This processing does not come without possible limitations. In particular, the finer resolution comes at the expense of (large) uncertainty in the estimation of the turbulence statistics.

Drones or unmanned aerial vehicles (UAVs) [120,121] offer a potentially powerful alternative to take airborne in situ measurements in complex terrain. Similar to full-scale airplanes, and with the same caveats, they can also be utilized to assess turbulence characteristics, e.g., [122]. Their size and flexibility would especially allow them to be flown close to complex and otherwise inaccessible terrain.

3.2. Air Quality Measurements

Operational air quality networks typically measure mandated air quality pollutants, such as nitrogen oxides, ozone, and particulate matter. Some selected sites also include some organic species and heavy metals. In Europe, about 8400 sites are continually reporting national air pollution data to the European Environmental Agency (EEA), which are archived in EEA data servers [123]. These datasets have been used extensively to validate air quality modeling systems, e.g., [124] and monitor the progress of the implementation of air quality standards. In mountainous areas, air quality stations are often located in urban areas along valleys. A limited number of mountain top observatories have provided valuable measurements of background conditions. The most prominent ones in the Alps include Sonnblick, Jungraujoch, Hohenpeissenberg, and Zugspitze observatories [125,126], while in the Apennines, Monte Cimone [127,128]. These observatories serve as European supersites and include an extensive suite of physical and chemical observations aimed at long term monitoring activities. They often provide the necessary space and logistics to conduct detailed atmospheric chemistry observations and deploy cutting edge technology, e.g., [129,130,131,132,133]. Long-term datasets at these observatories allow for the determination of trends in background air composition, against which models can be tested [134].

3.2.1. Past Projects

While it would go beyond the scope of this paper to provide an extensive review on community air quality campaign efforts in recent decades, we would like to highlight some significant undertakings in this area that were conducted in mountainous regions. The Megacity Initiative: Local and Global Research Observations (MILAGRO) field campaign was probably one of the most extensive air quality studies to date, conducted in Mexico City, which is situated in a basin 2240 m above sea level [135]. The basin is surrounded by mountain ridges on three sides, with a broad opening to the north and a narrower gap to the south-southwest, inducing characteristic valley and gap wind systems. The MILAGRO campaign included over 450 scientists from 150 institutions in 30 countries. The front range area in Colorado also has many decades of history of air chemistry campaigns. The area is impacted by complex flow patterns and gap winds. Poor air quality and associated airflows have been a major issue in Denver since the 1980s [136]. The most recent observational campaign investigating the impacts of the Denver Cyclone on regional air quality was the Front Range Air Pollution and Photochemistry Experiment (FRAPPÉ) [137,138]. The interplay between biogenic and anthropogenic pollutants in the Rocky Mountains was also a recent focus of Bio–hydro–atmosphere interactions of Energy, Aerosols, Carbon, H₂O, Organics and Nitrogen (BEACHON) [139,140]. In the European Alps, major air quality intensive projects include ALPNAP [10,11,12], MAP [84], Pianura Padana Produzione di Ozono (PIPAPO) [141], and Vertical Ozone Transports in the Alps (VOTALP) [19], while projects in the Pyrenees include Pollution Atmospheric in the Pyrenees (PAP) [142] and Pic 2005 [143]. These intensive observational programs provide the opportunity to study the complex processes of atmospheric chemistry and dispersion in much greater depth (e.g., ozone production efficiencies, organic chemistry, secondary aerosol production) than by using standard air quality observations. They also supply the necessary database to test atmospheric models performing intensive field measurements with multiple platforms not commonly adopted in traditional air quality networks.

3.2.2. Innovative Measurements

A new frontier in atmospheric chemistry observations are direct measurements of atmospheric fluxes based on the eddy covariance method. To date, measurements of organic and inorganic species in the gas phase, e.g., [144,145,146], as well as aerosols, e.g., [147,148], have been combined with micrometeorological methodology to study surface–atmosphere exchanges. These measurement setups have been deployed on the ground and on aircrafts in order to map biogenic and anthropogenic air quality pollutant and precursor emissions, e.g., [130,149,150,151,152,153,154,155,156,157,158,159].

Recent progress in low cost sensors (LCS), e.g., [160,161] presents new opportunities for air quality monitoring. Compared to traditional measurement systems, LCS devices tend to be less sensitive, less precise and less chemically specific to the compound or quantity of interest, but this tradeoff could be off-set by a potentially large increase in the spatial density of measurements that can be achieved by a network of sensors [162]. Currently, electrochemical sensors for NO_x and laser scattering sensors for particulate matter show some promise to become tractable monitoring devices in urban areas, typically exhibiting high concentration levels of pollutants [163]. Continuing progress in the development of LCS might lead to sensors that can also be deployed at rural sites. Pollution levels along Alpine valleys (especially the European Brenner transit corridor) are often high enough that even existing technology could potentially be integrated as part of a sensor network. Due to challenges in reproducibility and accuracy, LCS networks would currently need to be integrated within traditional air quality monitoring stations that can provide accurate reference measurements and help in discerning the cross-sensitivities of LCSs.

3.2.3. Remote Sensing

New opportunities to gain more comprehensive information of atmospheric pollutant distribution over complex terrain are also given by the increasing availability of different types of remote sensing techniques. For example, Raman and high spectral resolution (HSR) lidars are capable of retrieving profiles of aerosol backscatter and extinction coefficients. In addition, depolarization measurements are now possible, e.g., [164]. The 3 + 2 Raman lidar Polly [164], for example, combines three wavelengths and is able to obtain three backscatter and two extinction coefficients. The ability to observe depolarization allows for distinguishing certain aerosol types (e.g., mineral dust vs. smoke aerosol). Additionally, airborne aerosol lidars are a powerful means to map the spatial domain of fine–coarse aerosols (e.g., > 100 nm). Harnish et al. [165] analyzed vertical aerosol profiles in the Inn Valley from the airborne TropOLEX lidar [166], highlighting the complexity of the dispersion processes of pollutants in mountainous terrain. In particular, these aerosol lidar measurements provided an unprecedented three-dimensional picture of the distribution of air constituents in an Alpine valley, finding high aerosol backscatter intensities along sun-exposed slopes during wintertime. Similarly, Chazette et al. [167] in the French Alps attributed the vertical dispersion of aerosol layers to increased or decreased convection, depending on snow cover. In these examples, the combination of airborne lidar and in situ measurements proved to be a powerful solution to track pollution dispersion in complex topography.

Another approach to investigate the three-dimensional structure of air pollutants in complex topography is based on a combination of automated lidar ceilometers (ALC), sun photometers, and in situ aerosol measurements. Diemoz et al. [9], for example, used this setup to observe the diurnal and vertical structure of aerosol pollution transport from the Po Plain into the Aosta Valley. During pollution events, the ALC scattering ratio reached values > 30 with in situ PM₁₀ concentrations exhibiting values up to 100 μg/m³. Another possibility to obtain spatially resolved measurements of selected pollutants, such as NO₂, can be based on mobile Differential Optical Absorption Spectroscopy (DOAS), e.g., [168,169].

4. Meteorological Modeling for Pollutants Dispersion Applications over Complex Terrain

For adequately simulating both the advection and diffusion of pollutants, suitable information concerning both the mean flow and the turbulence statistics are needed for the dispersion modeling. These meteorological fields may be obtained either from measurements or from the output of diagnostic or prognostic models, or from a combination of both. However, in most cases, meteorological observations alone cannot provide the necessary details about the spatial variability of atmospheric variables, especially concerning turbulent quantities over complex terrain [170]. Similarly, when diagnostic models are fed by local observations only, they generally cannot produce good quality representative gridded meteorological fields in mountainous regions [171]. More reliable information about the mean meteorological fields can be obtained by using multi-sensor diagnostically derived nowcasting fields, such as the Integrated Nowcasting through Comprehensive Analysis (INCA) system [172]. However, the use of prognostic atmospheric models appears to be crucial for providing adequate meteorological input to dispersion models. Still, observations do represent a valuable integration into numerical prognostic simulations, through different data assimilation techniques, both as initial values and during the simulation. An in depth overview of the challenges related to data assimilation in complex terrain is provided by [173] in this Special Issue.

4.1. Meteorological Modeling Resolution

The requirements that should be satisfied by the output of prognostic models include many aspects. First, the meteorological model should be run at an adequate spatial resolution, in order to resolve at least the most important topographic features. In very complex terrain, this often implies performing meteorological simulations at sub-kilometer resolutions [14,174]. Indeed, it is well known that the model’s effective resolution is about six to eight times as large as the horizontal grid spacing [175]. The continuously increasing availability of computational resources allows the running of meteorological models at these resolutions, especially for the modeling of episodes [36]. Additionally, operational local area meteorological forecasts are nowadays performed with a horizontal grid spacing of the order of 1 km. This allows not only to resolve at least the most significant details in the orographic structure, but also to explicitly take into account an increasingly higher number of physical processes, generally leading to a more realistic reconstruction of the meteorological fields [176,177,178,179]. However, the high resolutions that can be now be achieved pose several other issues regarding the correct simulation of boundary layer processes over mountainous terrain. First, these resolutions often fall in the so-called gray zone or “terra incognita” [180], i.e., those scales where some atmospheric processes are neither sub-grid nor fully resolved, but rather partially resolved. This aspect is particularly important for turbulence parameterizations, which are crucial in pollutant dispersion applications. A comprehensive review of the problems connected to meteorological simulations falling in the gray zone is provided in [181] in this Special Issue.

4.2. Turbulence Parameterizations

Physical parameterizations implemented in meteorological models may suffer from intrinsic weaknesses when applied over complex terrain. In air pollution applications, this is particularly significant for the parameterization of processes occurring in the lowest atmospheric layers, where pollutants are mostly emitted and dispersed, such as turbulence and land–atmosphere exchange parameterizations. Turbulence parameterizations were derived using assumptions, such as horizontal homogeneity, which allows the treatment of turbulence as a one-dimensional vertical process, accounting for gradients only in the vertical direction, thus neglecting horizontal turbulence exchanges [182]. However, this is valid over flat terrain only; three-dimensional effects due to the presence of the orography in mountainous regions significantly affect turbulence statistics [183]. Some attempts have been made to extend turbulence parameterizations for applications over complex terrain [184], but further developments are still needed. Moreover, turbulence closures involve a series of coefficients, which are generally calibrated over flat terrain, and thus may not be appropriate for complex terrain situations [185,186]. These aspects can be of minor importance as far as only the mean fields over large scales and in simplified conditions (strong wind, weak convection, homogeneity), without any information about turbulence, are concerned, but become crucial in complex terrain, especially when meteorological simulations are used to drive dispersion models, where smaller scale motions, including turbulence quantities, are essential. To this regard, for air pollution applications, high-order turbulence closures should be adopted, in order to calculate the turbulence quantities needed for the dispersion model in a more physically based way [14,187]. Even if more sophisticated approaches, including higher order moments [187,188,189,190] and calculating turbulence transport [191,192], would be desirable, among the available options in the up-to-date numerical weather prediction (NWP) models, it is preferable to choose a model solving a dynamical equation for the TKE, since it provides more direct physical information [193]. This, however, implies considerably higher computational costs. Suitable turbulence information for the dispersion model could be supplied also by LES, which explicitly resolve the largest and most energetic scales of turbulence, while parameterizing the smaller scales, which are expected to be locally homogeneous and isotropic [194]. However, in many real cases, it is still not affordable to reach the spatial resolutions needed to resolve the most energetic eddies, the latter motions belonging to the upper part of the inertial subrange. In particular, under stably stratified conditions, which are among the most critical for air pollution dispersion, a resolution of few meters should be used [194].

4.3. Land Surface Parameterizations

Simulations of MoBL processes are also strongly affected by the performance of land surface models (LSMs), which calculate land–atmosphere exchanges of heat, moisture, and momentum. This is particularly true in complex terrain, where local circulations, controlling pollutant transport and diffusion are predominantly driven by the spatial heterogeneity of surface properties and surface fluxes. The performance of LSMs is determined both by the physical parameterizations adopted and by the quality of the static datasets, as well as of the initial conditions supplied. In fact, high-resolution simulations need high-resolution static datasets and initial conditions to take full advantage of the details potentially resolved by the model. While the available topographic datasets generally present suitable resolutions, the accuracy of other static datasets, affecting the behavior of LSMs—such as land cover, vegetation (including vegetation state), and soil properties—is generally lower. Moreover, land cover maps may be outdated, or the physical properties used for the characterization of the discrete number of land use categories may be not fully representative of real surface conditions. In particular, the representation of the land cover has been found to be an important potential source of error in the simulation of the surface energy budget, dynamically influencing the development of thermally driven winds and the evolution of cold pools and temperature inversions [195,196,197]. One of the most important parameters influencing the performance of LSMs is soil moisture, strongly affecting the partitioning of the available energy between sensible and latent heat flux, and thus the timing and strength of thermally driven circulations, as well as the MoBL height [174,198,199]. In fact, thermally driven circulations tend to be more intense, and boundary layer heights become deeper, over drier surfaces. However, initial fields of soil moisture (and of other soil fields) for mesoscale simulations are generally derived from the coarser grid of global models, which cannot represent the spatial variability typically occurring in complex terrain. Initialization problems can be overcome by coupling the meteorological model with a hydrological model [174], using the assimilation of data from in situ or remote sensing measurements [200], or performing a spin-up of the soil variables for an adequate time period. Spin-up can be performed by means of an offline land data assimilation system (e.g., the High-Resolution Land Data Assimilation System (HRLDAS), [201]), or by running the atmospheric model itself to cycle soil variables, e.g., [202,203]. Several studies have shown that spin-up can provide a realistic initial spatial distribution of soil moisture (and of other soil fields), reproducing more appropriately the variability of soil fields over complex terrain. This can affect the simulation of processes critical for the evaluation of pollutant transport and diffusion considerably, such as thermally driven circulations [204], boundary layer development and depth [205], and the near-surface temperature field [206], improving model results with respect to simulations with soil data simply interpolated from a coarser analysis.

A significant limitation of LSMs currently implemented in meteorological models is the fact that they are one dimensional, neglecting surface and subsurface lateral transport associated with topographic features [204]. This can be a relevant source of error where sloping surfaces are present, and, as a consequence, larger horizontal soil moisture redistribution is expected [207].

It is also known that LSMs often poorly simulate snow cover and snow cover/atmosphere exchanges [195,208]. An incorrect representation of snow cover may produce strong errors in the simulation of surface radiation budgets and consequently of near-surface temperature, thermally driven circulations, and the boundary layer height [209]. Since snow cover is often associated with persistent thermal inversions and cold pools, leading to critical air quality situations, this aspect is particularly important in air pollution applications. Similarly to soil fields, in mesoscale meteorological simulations, snow cover is usually initialized with interpolated data from global models, suffering from analogous weaknesses as those illustrated for soil properties, as snow cover over mountainous terrain presents a high spatial variability, strictly connected to altitude variations. Moreover, reliable observations of snow cover and depth for data assimilation are scarce. However, satellite products, such as those provided by MODIS (Moderate Resolution Imaging Spectroradiometer)/Terra [210] can now be used to improve the snow representation in numerical models [36].

4.4. Urban Areas

Meteorological models should also be able to simulate the microclimatic modifications induced by urban areas, and in particular the local alterations of the wind field and of the thermal stratification, to correctly reproduce the dispersion of pollutants emitted in cities. In recent decades, specific schemes have been developed to parameterize surface exchange between the urban environment and the atmosphere and have been implemented in meteorological models [211]. Several studies showed that their application considerably improves the simulation of meteorological fields in urban areas, especially when detailed and high-resolution urban morphology datasets are available [212].

5. Pollutant Dispersion Modeling over Complex Terrain

5.1. Input Data for Dispersion Models

Dispersion modeling in complex terrain requires, in principle, the same meteorological input as anywhere else in the world (see Section 3.1). However, as highlighted in Section 4, a critical challenge in complex terrain arises from the fact that the available local observations may not be representative at all of the dispersion processes to be modelled, and certainly, spatially distributed input is required. Indeed, even when using advanced input data, the meteorological information passed to the dispersion model may contain significant errors, introducing inaccuracies both in mean and turbulent quantities. For example, Rotach et al. [91] have compared wind direction from different approaches (prognostic, with a grid spacing of 1 km; or diagnostic, down to a 50 m grid spacing) with a number of relatively close (within a few kilometers) observational sites in the Inn Valley (Austria) and found root mean square differences of more than 45° (for deliberately chosen challenging situations like foehn or a frontal passage). Goger et al. [213] have used data from the same area to investigate the uncertainty in the simulation of TKE using a numerical model with a 1 km horizontal grid spacing, and reported a root mean square error of a few tenths of a m² s⁻² (0.36–0.51, depending on the closure details in the model) for several selected days. A possible alternative, in view of quantifying at least this uncertainty, could be represented by ensemble dispersion modeling approaches. However, such an approach has not yet been tested to the knowledge of the authors, or at least, not in complex terrain.

Focusing on the input data concerning turbulence, most dispersion models will need velocity statistics (such as velocity variances, or possibly higher order moments like skewness, e.g., [214]). The type of dispersion model and its turbulence closure may make it necessary to add one or the other variable, such as the dissipation rate of TKE. Turbulence variables are commonly obtained by using a meteorological pre-processor (i.e., similarity relations applied to mean meteorological data to retrieve the turbulence statistics), based on data from a nearby weather station or from a numerical model, with the same problems of representativeness as for mean quantities. Moreover, the calculation of dispersion parameters required for dispersion models is generally performed by means of oversimplified parameterizations, again developed for flat terrain applications, adding further uncertainty to the dispersion simulations. In particular, dispersion models generally use dispersion parameters derived from scale analyses of the NWP model surface layer, e.g., [215], provided by the NWP model, or internally calculated from the mean fields in the dispersion model. However, if some turbulent information is available from the NWP model, an alternative approach can be adopted, accounting for these variables. In this way, more physics is introduced in the diffusion coefficient, with, in principle, more reliable results [14].

5.2. Eulerian vs. Lagrangian Models

The two main alternative approaches in the numerical modeling of atmospheric pollutant dispersion are the Eulerian approach, where the frame of reference is fixed with respect to the Earth, and the Lagrangian approach, using a coordinate system moving with the average atmospheric motion. Both Eulerian and Lagrangian dispersion models are three-dimensional (3D) prognostic tools, taking into account the variability in meteorological fields and the effects of inhomogeneities due to terrain complexity, land use heterogeneity, and even the presence of obstacles [216]. In the grid-cell Eulerian models, the 3D advection–dispersion equations are numerically solved on each grid point at each time step, and chemical reactions can be resolved. Therefore, the concentration of the substances of interest is advantageously estimated together with the meteorological variables. Eulerian models are effective when the gradients of concentrations can be well resolved by the numerical grid. This implies that, while the pollutant dispersion is well described when the cloud size grows to the scale of the grid cell, close to the emission point and in the first phase of the dispersion process their ability is affected by the limit of the grid resolution, being unable to capture the sub-grid diffusion. Similarly, the grid size should be fine enough to adequately resolve the interaction between the plume and the orography. Depending on the scale of interest and the orography complexity, very high grid resolutions may need to be adopted for properly detailing the dispersion scenario, with an increase in the computational cost.

On the other hand, the Lagrangian approach is grid-free, it follows, at all scales, the motion of individual plume parcels, whose paths are modeled based on a random walk process, and it describes all phases of dispersion with the same accuracy, including the near-source region. In Lagrangian particle models, the dispersion of the airborne pollutant is simulated through fictitious particles, each containing a small amount of the emitted tracer mass. The particle size is small enough to move according to the smallest eddies and large enough not to be influenced by the viscosity. The local wind drives their mean motion and the diffusion is determined by velocities obtained as the solution of Lagrangian stochastic differential equations, providing the statistical characteristics of turbulent flows. Different portions of the emitted plumes can experience different atmospheric conditions, allowing a realistic reproduction of the complex atmospheric phenomena even in mountainous regions. The theoretical framework on which modern Lagrangian models are based derives from the seminal work of Thomson [217]. He established a fundamental criterion to be met by Lagrangian models, the well-mixed condition, stating that, with reference to the position–velocity phase space, particles initially well-mixed must remain so for all subsequent times. Lagrangian models are pure dispersion models and require meteorological variables that are generally provided as 3D-gridded fields and produced by diagnostic or prognostic Eulerian models. Thus, the resolution and the accuracy of the input fields condition the level of detail in the simulation of the dispersion processes. A mesh is generally also used for calculating the concentration in the domain of interest, by computing the statistics of the particle trajectories. Advanced Lagrangian particle models, when including both vertical and horizontal derivatives in the stochastic equations for the wind velocity, could benefit from the meteorological fields produced by Eulerian LES models. However, in order to treat the 3D non-Gaussian velocity statistics needed in the MoBL, skewed non-Gaussian turbulence has to be solved by the Lagrangian model. The Thomson [217] well-mixed condition provides a unique solution in one dimension and thus a unique solution for the 3D case would require that variables in one dimension are independent from those in the other dimensions. Inhomogeneous non-Gaussian turbulence was treated in 1D [218] and 2D models [219]. Approaches to obtain a 3D skewed turbulence have been proposed (see, for instance, [220,221,222,223]), yet they are not generally operational.

First [224] and second order [225,226,227,228,229,230] chemical reactions can be solved in Lagrangian models, but only a limited number of chemical equations can be taken into account. It is generally recognized that the segregation of the chemical reactants cannot be neglected in short-term concentration predictions [231] when chemical reactions take place before the pollutants are well mixed by the turbulence. Therefore, a correct description of this phenomenon needs the estimation of the concentration fluctuations. In particular, when the ratio of the chemical reaction rate to the turbulent mixing rate, the Damköhler number, is small and the equilibrium is almost instantaneously reached, chemistry is not affected by the turbulent dispersion, but when the Damköhler number is ∼1 or larger, the segregation of the reactants may be important [231]. In this case, concentration fluctuations are needed and Lagrangian models able to predict concentration moments of orders higher than the first (see, for instance, [228,232,233,234,235,236,237]) can be successfully used. On the other hand, Eulerian models provide only mean concentrations and cannot take into account this effect. Regarding complex terrain, where turbulence scales are generally smaller, the Damköhler number is likely to be higher, and hence the segregation of the chemical reactants should be considered for longer times. In other words, the time range in which segregation cannot be neglected becomes longer when reacting plume dispersion takes place in complex terrain.

5.3. Small-Scale Features Matter

As highlighted in Section 2, complex terrain affects the transport and dispersion of pollutants through both mechanical and thermal effects. Obstacles, like hills or mountains, modify the flow and hence the motion of the plume emitted from a source. For example, depending on the horizontal scale of the obstacle, the plume can either climb over or turn around the relief. All these conditions represent a challenge in modeling pollutant dispersion; the traditional codes, developed for flat terrain, cannot be applied, but new approaches should be investigated [14,186,238,239,240,241]. In particular, models adopted for regulatory purposes may fail in correctly capturing pollutant dispersion in complex terrain. Tomasi et al. [14] highlighted that the Gaussian puff model CALPUFF, widely used for regulatory modeling applications, poorly represented ground-level concentrations of a tracer released by a chimney in the Adige Valley (Italy) and in particular its interaction with the orography, whereas better results were achieved with the Lagrangian model SPRAY-WEB, especially when TKE information from the meteorological model was used (Figure 4).

5.4. Low Wind Speeds and Meandering

A challenge for pollutant dispersion models is determined by the proper description of meandering; in these situations, the problem of identifying and defining a precise mean wind direction makes the simulation of dispersion rather difficult. Indeed, high air pollution concentrations are often associated with low wind speeds, but due to meandering, the ground level concentration can be rather lower than what is predicted by numerical models not accounting for, or poorly resolving, these particular features. In Lagrangian stochastic dispersion models, different solutions have been proposed, for example: an “ad hoc” algorithm based on a fluctuating plume model [242], a Markov chain Monte Carlo model simulating negative lobes in the Lagrangian autocorrelation function [243], a system of two coupled Langevin equations for the two horizontal wind components explicitly accounting for meandering [244], and an analytical correction factor based on mass balance adjusting the modeled point source concentrations at the lowest model level for concentration gradients [245]. A more detailed discussion on the methods used in Lagrangian models for light winds can be found in [246].

5.5. Plume Rise

One of the key factors for the correct estimation of the dispersion of airborne pollutants and for the evaluation of ground-level concentrations is the computation of plume rise. In the current operational models, e.g., [247], the plume rise is computed assuming an air parcel’s rise based only on the buoyancy terms [248], using the heat release, the wind velocity, and the friction velocity during the day and the static stability at night. The plume is injected at the final plume height from the center of each grid cell that includes the source location. These assumptions can lead to poor approximations, particularly in complex terrain, where the plume does not move horizontally but follows the topography. Furthermore, plumes with a large buoyancy flux are likely to reach the top of the boundary layer during the day and to partially penetrate above the capping temperature inversion layer [249]. In the Eulerian dispersion models, the calculation of plume rise is based on the fluid dynamic equations, namely on the mass, momentum, and energy conservation equations [250]. In Lagrangian particle models, plume rise can be dynamically computed, i.e., each particle, at each time step, can respond to local conditions, such as wind speed and direction, ambient stability, and turbulence, e.g., [251,252,253]. This makes Lagrangian particle models more suitable for correctly simulating the plume rise in complex terrain, where the assumption of horizontal homogeneity cannot be made.

5.6. Odor Dispersion

A particular problem related to atmospheric dispersion is the assessment of odorous substance transport. In several countries, odor assessment is based on the so-called odor hours, defined by at least 6 min of perceivable odor concentrations in an hour. Modeling odor hours requires the determination of the 90th percentile of the corresponding cumulative frequency distribution. It can be determined by assuming a probability density function (PDF) that is generally completely defined by a concentration fluctuation intensity i and by a mean concentration [254,255]. Obtaining an estimation of i requires a more sophisticated approach. Many models have been developed for the second and higher order moments of the concentration PDF over the last ten years. Thomson [256] and Borgas and Sawford [257] proposed the two-particle model, Pope [258] and Cassiani et al. [259] proposed the PDF method, while the fluctuating plume model was suggested by Ferrero et al. [230], Luhar et al. [232], Franzese [233], Mortarini et al. [234], Bisignano et al. [235], Yee et al. [260], and Yee and Wilson [261]. While the two-particle model can be applied only in homogenous and isotropic conditions, the fluctuating plume model is able to take into account the flow inhomogeneity and, as a consequence, is suitable for applications in complex terrain. Recently, Manor [236] used a Lagrangian particle model to simulate the concentration variance dispersion and Ferrero et al. [237] suggested a formulation for the variance dissipation time scale to be used in the same model. This approach was used by Oettl and Ferrero [254] and Ferrero and Oettl [262] to assess odor dispersion. For a complete review on this subject, see Ferrero et al. [263].

6. Conclusions

The paper reviewed open questions and challenges posed by both the measurement and modeling of processes controlling the advection and diffusion of air pollution, with a particular focus on the additional challenges occurring in complex terrain.

Our capability of monitoring and predicting air quality in mountainous areas still suffers from significant limitations. Some of them are inherent to the intrinsic complexity of atmospheric processes over the variety of landforms, and some depend on the limitations that impede the adoption of conventional measurement systems or conventional modeling schemes originally developed for simpler terrains. Therefore, various challenging questions still await adequate progress, not only for researchers in the field, but also, and even more, for authorities in charge of managing air quality in mountainous regions, relying on appropriate methods.

As outlined in the paper, technology nowadays offers more and more sophisticated instruments and for progressively decreasing costs. Given the large spatial variability of meteorological variables over complex terrain and the logistical difficulties of operating dense networks of surface weather stations, an important advancement for monitoring both mean and turbulent meteorological quantities will be represented by the increasing accessibility of remote sensing instruments and the related development of suitable methods for the interpretation of the measurements. The availability of spatially distributed meteorological fields is in fact crucial for both the initialization and the validation of meteorological models. In particular, spatially distributed information on turbulent quantities would be highly beneficial both for the understanding of atmospheric dispersion processes and for the evaluation of the parameterizations implemented in meteorological and air quality models. However, suitable techniques and methods to derive turbulent variables from remote sensing measurements in complex terrain are, in many cases, still not consolidated and exhibit specific challenges and limitations. Therefore, further efforts in this area of research would be highly desirable.

Similarly to meteorological fields, dense networks of air quality stations should be implemented over complex terrain to capture the spatial variability of pollutant concentrations, to detect possible critical situations, and to provide observational evidence for the evaluation of dispersion models. Operational air quality networks usually do not satisfy these requirements and, as a consequence, reliable datasets to test dispersion models over complex terrain are scarce. Suitable tracer experiment strategies should then be promoted to support model validation efforts and to provide datasets to modelers that can be used as benchmarks for testing dispersion models in complex terrain [13]. A promising way to establish denser permanent air quality networks, or to complement intensive field measurements, is the progress in the accuracy of low-cost instruments, even if challenges in reproducibility and cross-sensitivity still hinder their widespread use. Traditional air quality measurements can also be complemented by remote sensing instruments, which are increasingly adopted to provide spatially distributed fields of air quality parameters. A large-scale observational effort should ideally combine all these measurement strategies to gain the most complete picture of air chemistry in mountainous environments. A valuable example in this context is the newly established Innsbruck Atmospheric Observatory (IAO, [159]), which aims to provide basic research infrastructure, along with cutting edge atmospheric chemistry and physics observations to obtain an in depth view on air quality constituents and transport in the Inn Valley (Figure 5). Such comprehensive datasets are, however, still often missing in complex terrain, which significantly impacts the predictability of air quality across mountain ranges.

As for modeling applications, we have highlighted that, in complex terrain, the different steps involved in the simulation of pollutant dispersion—the creation of the meteorological input and its implementation in the dispersion model, dispersion modeling—are affected by additional challenges and uncertainty when compared to simple terrain. We have argued that sub-kilometer meteorological simulations are often needed over complex terrain to correctly capture at least the mean meteorological fields and thus to provide an adequate input into dispersion models. However, the applicability of state-of-the-art numerical weather prediction models, based on Reynolds averaging, is questionable at these resolutions, falling in the so-called “terra incognita”. Alternative solutions could be represented by large eddy simulations, but the applicability in real cases over complex terrain is still limited to specific situations and for research applications, due to their high computational cost. For operational purposes, a resolution of the order of 1 km is nowadays commonly adopted and this will probably remain the standard in the coming years. It follows that advancement in the parameterization of turbulence processes for mesoscale model simulations is still needed, especially over complex terrain, where most of the hypotheses commonly made for their development often fail. In particular, efforts in the representation of horizontal inhomogeneity in planetary boundary layer parameterizations should be pursued, implementing 3D turbulence schemes with different horizontal and vertical length scales accounting for turbulence anisotropy [183,213]. Moreover, high-resolution simulations in complex terrain often suffer from inadequate initial conditions regarding soil parameters and snow cover, which deeply affect near-surface meteorological fields and the development of the boundary layer. Suitable techniques for initializing these variables should then be developed and adopted to take full advantage of the details potentially captured by high-resolution simulations. In some situations, resolutions higher than what can be achieved by mesoscale models, especially for operational needs, would be preferable to better resolve the interaction between the flow field and the orography. For this purpose, diagnostic flow models, specifically developed for complex terrain applications and relying on suitable data assimilation techniques, might be a computationally cheaper solution to provide high-resolution mean meteorological fields to dispersion models. The input needed by dispersion models includes, not only mean variables, but also turbulent quantities. These are commonly obtained by using a meteorological pre-processor, which can ingest data from observations or from a meteorological model. When suitable turbulence information is not available in the input data, the calculation of dispersion parameters is performed from the mean meteorological input by means of parameterizations developed in the framework of MOST, which again is not strictly applicable in complex terrain. This introduces further approximations, which can be at least partly reduced if some data about turbulence (for example, TKE) are produced by the meteorological model by using higher order turbulence parameterizations.

Similarly to meteorological models, dispersion models should be able to account for the enhanced spatial variability typical of mountainous regions and, in particular, for small scale turbulence. It results that Eulerian models should be applied at very high resolutions in complex terrain to adequately capture the interaction between pollutants and the orography, with the same computational limitations of meteorological models. As a consequence, Lagrangian particle models represent a more suitable tool in these conditions, at least for research applications. Some Lagrangian solutions have been proposed in recent years to simulate simple second order chemistry, but more research is needed in this particular field to extend the applicability of Lagrangian models to a more complex reactive pollutant chemistry. However, differently from Eulerian models, Lagrangian models are able to account for the segregation of the reacting pollutant, which may be more remarkable in complex terrain with respect to flat terrain.

This review highlighted that, despite past projects and research efforts greatly contributing to the progress of our knowledge in this field, several aspects related to the understanding, observation, and simulation of dispersion processes over complex terrain still remain to be investigated. Therefore, further research activities should be pursued in the directions outlined in this paper, including large community projects, in order to take full advantage of the opportunities offered by the continuous technological advancement of both observational and computational infrastructures.