Impact of Saharan Dust Intrusions on Atmospheric Boundary Layer Height over Madrid

1. Introduction

Aerosols constitute a public health risk, as they can be considered one of the main factors contributing to poor air quality [1]. Aerosols are mainly concentrated within the atmospheric boundary layer (ABL), with long-range transport layer sometimes located in the free troposphere. This heterogeneous distribution make their characterization difficult The ABL is defined as the lowest part of the atmosphere, influenced by the Earth’s surface by means of exchange of energy and moisture [2]. It plays a critical role in air quality forecasts [3] and greenhouse gas concentration budgets [4] and it is mainly characterized by turbulent processes. The influence of solar radiation creates a daily evolution cycle of the layer [5]. The cycle, in clear-sky situations, starts with the increase of ground surface temperature after sunrise, which intensifies the convection. It produces ascension of warm air masses and downward displacement of colder air masses, which creates a growing mixing layer (ML) [6]. It is called this way because substances emitted into this layer disperse gradually horizontally and vertically due to the turbulence. When sufficient time is given and there are no sinks, the ML become completely mixed [4]. The vertical turbulent mixing processes produces a strong aerosol gradient between the ABL and the free troposphere. When the ML is well developed, usually at the end of the daytime period, the ABL height (ABLH) is estimated similar to the ML height (MLH), as the turbulent mixing produces a nearly homogeneous distribution of aerosols along the complete layer. Later in the day, the gradual reduction of incoming solar irradiance during the early evening transition period causes a weakening of the turbulence and the convective processes. This creates a nocturnal boundary layer (NBL) close to the surface, that is stable and stratified. The remnant of the daytime ML is located above the NBL and it is called the residual layer (RL). When the sun rise again, a new ML begins to grow rapidly, eroding firstly the NBL and then entraining into the RL, producing an early morning transition. Since different surfaces respond differently to the solar heating, the development of the ABL is influenced by the surface albedo of the underlying surface [7]. Therefore, the combined effects of the synoptic atmospheric conditions, such as atmospheric stability or wind shear, and surface characteristics such as cover, roughness or topography, determines the horizontal variations in ABL dynamics [8]. There are a wide range of applications with high societal, economic, and health impacts, such as air quality [9] the generation of renewable energy [10] or numerical weather prediction [11] that benefits from a better understanding of these ABL processes.

The traditional methods employed to characterize the ABLH, radiosounding in the frame of the World Meteorological Organization Radiosounding Global Network [12]; present spatial and temporal limitations as they are usually launched only twice per day in most airports worldwide. This scarce observations prevent an accurate determination of the ABLH variations, both temporal and spatially, compromising their representativeness for urban and regional scales. It was documented that the calculation of the MLH based on potential temperature profile or Richardson number performed by numerical meteorological models produces more than 50% uncertainty in shallow boundary layers and 20% in deeper boundary layers due to scarcity in the input data [12].

Continuous profiling of the entire ABL vertical extent is nowadays possible thanks to recent advances in ground-based remote-sensing technology and algorithm development [13]. The high temporal and vertical resolution of remote sensing instruments permits a precise automatic detection of ABL sub-layer heights [14]. Among novel remote sensing methods, a promising one is the ceilometer, a particular type of lidars operating with a single-wavelength and originally intended for cloud base height detection [15]. Currently, ceilometers provide continuous high resolution aerosol backscatter profiles (every 15 s) with good spatial resolution (tenth of meters) and a large vertical range (up to several km ) in unattended continuous operation. Furthermore, ceilometers are usually operated in networks, such as EUMETNET EPROFILE [16] and ICENET (Iberian Ceilometer Network, [17]).

The ABLH retrieval algorithm for these instruments is based on the aerosols vertical profile. It assumes that aerosol concentrations are lower in the free troposphere than in the mixing layer, producing a strong negative gradient clearly observable in the backscatter profiles [18]. The high temporal and spatial resolution allows the study of the aerosol concentration fluctuations produced by the constant interexchange of airmasses, those polluted with aerosols within the ML moving upward are exchanged with those clean moving downward from the free troposphere. A disagreement can result between the height of aerosol layers derived by lidar profiles and the radiosounding-derived MLH due to the inconsistency between the thermal profiles and the aerosol profile, especially during morning or evening transitions [19]. ABLH retrieval algorithms take advantage of these characteristics to determine the height, such as gradient method [20], wavelet covariance transform [21] or edge detection method [22]. More advanced methods extend the analysis to two dimensions (temporal and vertical) in order to guarantee temporal consistency [23], or applies graph theory to track the diurnal evolution [24]. A reliable new method called STRATfinder [25] combines these last two methodologies applying a backward propagating layer, from the end of the day, in order to determine the type of layer by minimizing a cost function from the forward and backward runs. This layer attribution is normally assisted by commonly available surface measurements of radiation and temperature in order to decrease its uncertainty [23].

However, more complex atmospheric situations can jeopardize the observable strong negative gradient in the backscatter profiles. For instance, desert dust intrusions, volcanic eruptions or forest fires [26] can produce long-range transport of aerosols that affects the vertical distribution over the site due to the presence of lofted layer, which can either mix with the ABL or remain above it for long periods. This advection of aerosols increases the layer retrieval uncertainty, challenging any aerosol-based method. [27]. The Iberian peninsula is located close to the Sahara desert and regularly receive desert dust events, occurring when the dust is injected and transported throughout the atmosphere over long distances. On a global scale, the apportion of desert dust have been estimated up to 40% of aerosol mass yearly injected into the troposphere [28] and in particular Sahara desert emit half of the world atmospheric mineral dust [29].

The main objective of this work is the assessment of the impact of Saharan dust intrusions on the ABLH using ceilometer signals, along four years period 2020-2023. The database of ABL heights of continuous measurements have been classified regarding the most frequent patterns of synoptic circulation. The weather conditions over the Iberian Peninsula were analyzed using cluster analysis of sea level pressure fields and six typical synoptic meteorological patterns (SMPs) were identified. Section 2 describes the methodology employed including the instrumentation (ceilometer profiles) and algorithms (STRATfinder), the synoptic meteorological patterns and the methodology followed to identify Saharan dust intrusions. Section 3 summarizes the main results obtained when the datasets are differentiated by season and by synoptic meteorological patterns and Section 4 discuss the main findings of the work.

2. Materials and Methods

2.1. Experimental Site

The CIEMAT (Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas, 40.457◦N, 3.726◦W, 669 m asl) site is located Madrid, the capital of Spain placed roughly at the center of the Iberian Peninsula (See Figure 1). The site may be considered as an urban background station since there is a park (dehesa de la villa) closeby and the main traffic avenues are located moderately far away. Regarding the geography of the site, there are mountain ranges to the north-northwest (Sierra de Guadarrama, with maximum altitude of 2420 m. agl.)and also to the south (Montes de Toledo). the Madrid city’ population reached nearly 6 million inhabitants considering the metropolitan area and surrounding small towns, being one of the most populated regions in Spain. Main contributors to the Madrid air pollution plume are traffic emissions and residential heating, so it can be considered as typically urban. There is no significant industrial activity, only light factories [30]. Other contributions to the Madrid pollution be reduced to long-range transport of dust from the Sahara desert. In certain meteorological situations, Saharan dust intrusions can significantly contribute to the aerosol concentrations measured in the Madrid region [31]. The arrival of Atlantic air masses produces a cleansing effect on the Madrid atmosphere, with a significant reduction in the particulate matter levels. The weather in Madrid is described as continental Mediterranean, with prevalent clear-sky conditions, hot dry summers and cold winters [32]. During a large part of the year, the atmospheric situation is governed by the Azores high-pressure system which can produce periods of stagnation with high stability, poor ventilation and increases in local air pollution in winter.

2.2. Instrumentation

At the site, a ceilometer Lufft CHM15k-Nimbus belong to the instrumentation employed in the MDR-CIEMAT ACTRIS station [33], and, since 2020, the ceilometer data is available at ICENET [17]. The instrument operates with a pulsed Nd: YAG laser with output power per pulse of 59.5 mW at 1064 nm and a repetition frequency ranging between 5 and 7 kHz. The overlap between the telescope (mounted in biaxial configuration) and the laser beam is 90% complete between 555 and 855 m above ground level, according to the instrument datasheet, but Heese et al. [34] found that the complete overlap is located at 1500 m a.g.l. To improve this characteristic, a correction function was introduced to reduce the incomplete overlap producing a useful signal down to 232 m agl, where the laser beam divergence is 0.3 mrad and the telescope field of view, 0.45 mrad [36]. This overlap function has been successfully tested in previous works [33,35]. The ceilometer files for complete day are processed by means of the STRATfinder algorithm, available under the GNU General Public License v3.0 [25]. This recently developed algorithm estimates an auxiliary layer height that is tracked backwards in time from midnight to noon, then the Dijkstra algorithm [37] is applied to track MLH, ABLH and the mentioned auxiliary layer providing as output the estimations of the MLH and ABLH. A fast Fourier transform function is applied to exclude non-reliable estimations of MLH and ABLH from attenuated backscatter profiles. This also filters cases of rain, snow and low clouds, so these situations are excluded from the dataset. As it was mentioned in the introduction, ABLH estimations are more challenging during the evening decay of the mixed layer, so the algorithm assist in the identification by connecting individual paths merging the auxiliary layer and the preliminary ABLH to provide the final MLH and ABLH estimate for the whole 24 h period. The atmospheric vertical profiles are obtained from the backscattered signals with a temporal resolution of 15 s and a spatial resolution of 15 m, and averaged to 1 hour before stored in the database as daily files. For each day, the MLH at 16:00 was selected as representative height for that day. This allows a more reliable comparison among days than selecting the maximum value for each day, as it was tested in the data processing. The final dataset comprises a total of 1461 days between January 2020 and December 2023.

2.3. Methodology for Classification of Days Based on Meteorological Fields

The most frequent synoptic meteorological patterns (SMPs) over the Iberian Peninsula were identified using sea level pressure (SLP) at 12 UTC over the period of 2001–2019 [38]. The non-hierarchical k-means cluster analysis was employed to assign each day to a specific SMP. The process consisted in 6 stages: Firstly, the reanalysis global fields of SLP for each day at 12 UTC were downloaded from the NCEP/NCAR (NOAA/OAR/ESRL PSD, Broadway Boulder, CO, USA) reanalysis dataset open access repository [40]. Then, the initial conditions were established selecting several k SLP fields, based on results concerning typical SMPs occurring in the western Mediterranean basin, which were obtained in previous studies. Thirdly, the Euclidean distance for the value of every grid point was calculated between each SLP data field and each k-cluster center was calculated and summed. The arithmetic mean of all their members, grid point by grid point provides new cluster centers, once every SLP data field was grouped into a specific cluster. Repeating this iteratively, a stable solution can be reached when each SLP field remains in its cluster from one iteration to the next. The results are composite maps representing the SMPs (see Figure 2) obtained by averaging all the SLP fields grouped into each cluster. The final number of clusters (k) was selected using variance plots. As a validation procedure, a complete dataset of parameters were analyzed for assuring the physical meaning of the SMPs. Measurements of meteorological variables performed in an instrumented tower, atmospheric stability parameters from numerical models and radiosondes were employed, which characterize the daily meteorological situation at the surface level in the central region of the Iberian Peninsula. This procedure provides physical meaning to the six different and realistic SMPs over the Iberian Peninsula. This classification have already been used in former studies to evaluate the influence of atmospheric conditions on airborne allergenic pollen and spores in Barcelona [40] and on the temporal evolution of the mixing layer over Madrid determined by different methods [41]. For this study, the classification can be applied to the period 2020-2023, assigning each day to a specific SMP for the subsequent analysis of the estimation of the ABLH,

2.4. Identification of Saharan Dust Events

The occurrence and duration of Saharan dust episodic days over the Iberian Peninsula was determined by a robust method based on four different components. Firstly the daily interpretation of air mass back trajectories computed by the HYSPLIT model [42], secondly, synoptic meteorological charts, thirdly, maps of aerosol index of Ozone Monitoring Instrument-OMI and NASA SeaWiFS images and finally, three dust forecast models outputs, namely: SKIRON-University of Athens (http://forecast.uoa.gr), MONARCH-Barcelona Supercomputing Center (https://dust.aemet.es/products/daily-dust-products) and NAAPS-Naval Research Laboratory (NRL), Monterey, CA (http://www.nrlmry.navy.mil/aerosol/). Once this episodic days are identified, the dust contribution to the PM10 daily records can be quantified performing an statistical analysis of the time series of PM10 values registered at regional background monitoring sites, classifying the intensity of the dust episodic day by their dust contribution levels. Different approaches has demonstrated the feasibility of this method [43] and nowadays is the Spanish and Portuguese reference method to identify and quantify African dust contributions to PM10 levels since 2004. Currently, this is one of the official methods recommended by the European Commission for evaluating the occurrence of African dust intrusions and quantifying its contributions [44]. In this work the African dust events were also separated by dust load (concentration of dust in the PM10 fraction registered at the rural background station “El Atazar” (Air Quality Network of the Madrid Region: https://gestiona.comunidad.madrid/azul_internet/run/j/AvisosAccion.icm during Saharan dust episodic days) considering clean days those without contribution of dust from the Sahara, low dust days those with concentration below the median plus half the standard deviation of the series, that yields a value of 30 µg/m³, and the high dust days are those above that concentration.

3. Results

3.1. Daily Evolution of the ML, Quicklooks and STRATfinder Estimations

Figure 3 shows two cases as example of the STRATfinder estimations. Figure 1a shows a usual daily evolution case, occurred in 21 June 2022. The MLH (black circles) shows low values during the night and part of the morning, close to the minimum (232 m agl). During the same period, the ABLH (red crosses) reaches 1.8 km and is detected at the layer edge between the ABL and the free troposphere, indicating the residual layer from the previous day. The convection starts at sunrise (5:59 am), observable in the figure as vertical cyan lines along the day, detecting the strong upward movement of airmasses, up to 3 km at 18:00, just before sunset. Both estimations, MLH and ABLH are equal at 10:00, indicating a complete growth of the ML over the RL, although later on, the algorithm estimates lower MLH from 11:00 to 15:00, probably due to some inner structure in the aerosol layer. The agreement remains between 15:00 and 17:48, sunset time. After that time, the algorithm identify the MLH as the first detected layer close to the ground, meanwhile the ABLH tracks the residual layer that slowly transitions into the next day. Such daily evolution occurs consistently on clear-sky days, especially in summer. As it was mentioned in section 2.2, the MLH at 16:00 will be selected for this day, in this case, close to the maximum and in coincidence with the ABLH.

Figure 3b shows a case with a strong Saharan dust intrusion arriving at the site. As it can be seen, the daily evolution is different, with a large ABL, reaching up to 4 km, and a ML showing a delayed development during the morning, caused by the effect of the dust-rich airmass over the radiation that reach the surface and also over the meteorological conditions. Aloft aerosol-rich layers are observed between 9:00 and 21:00 from 2 to 3 km , with clouds at the end of the day located at 4 km, on top of the ABL. For this day, the selected MLH at 16:00 reached lower value than the ABLH due to the effect of the aloft dust-rich layer. The study of this influence of the presence of Saharan dust on the MLH is the main aim of this work, so the database have been discriminated between days with Saharan dust influence and clean days.

3.2. ML Heights Estimation for 2020– 2023

Figure 4 shows all the data, both MLH (black crosses, left axis) and daily dust loads (orange circles, right axis) for the four years studied (2020 - 2023). Several features can be observed, namely a clear seasonal behavior of the MLH along the four years, with maxima in summer, around 4 km, and minima in winter, less than 300 m and with a large variability. Development of Saharan dust episodes show the same seasonal behavior, with the highest number of episodic days in summer (40% of the total) and the lowest in winter (16% of the total). Finally, 361 days of the period 2020-2023 were categorized as Saharan dust episodic days (25% of the total). On average, daily dust load reach 18.1 µg/m³ along the study period with the highest mean values registered in winter (24.4 µg/m³) and spring (21.2 µg/m³). This is mainly a consequence of the occurrence of intense episodes with very high dust load, higher than 100 µg/m³ and even reaching 340 µg/m³ on Feb 2022, occurring in February and March. Considering the threshold value of 30 µg/m³ of daily dust load, the database identify 50 high and 311 low Saharan dust episodic days. Such discrimination will allow a more detailed study of the influence of Saharan dust, as it will be explained below.

A first analysis of the complete database (all days considered) shown in Figure 4 was attempted separating by season. The results, plotted as boxplot graph (Figure 5), show median MLH values of 721 ± 502, 1694 ± 628, 2262 ± 875 and 1371 ± 638 m for winter (DJF), spring (MAM), summer (JJA) and autumn (SON) respectively, quantifying the seasonal effect mentioned before. The season are selected using complete months, corresponding winter to December, January and February, or DJF, spring to March, April and May, or MAM, summer to JJA and autumn to SON:

More detailed analysis can be accomplished by classifying days by the synoptic meteorological patterns, as explained in the methodology section 2.3. Figure 5b shows that the highest MLH corresponds to SMP2 (2004 ± 866 m), that mainly occurs in summer (53% of cases), while the lowest MLH correspond to SMP1 (630 ± 460 m) a clear winter type (86%), without cases in summer. The other synoptic patterns show a similar trend, with the next lowest cases corresponding to SMP6 (1002 ± 593 m), the other main pattern for winter (50%), while the other three show increasing MLH with the frequency of their cases in summer (SMP3: 980 ± 677 m (3% of cases in summer), SMP4: 1541 ± 721 m, 16% and SMP5: 1802 ± 804 m, 32%).

In order to characterize the different SMP, the distribution of cases for each season respect to the SMP are represented in Figure 6. The number of cases are converted into frequency by dividing them for each season by the total number of cases assigned to that SMP. As it can be seen in the figure, a clear representation of SMP1 cases in winter (86%) and none in summer, and a strong representation of cases of SMP2 in summer (53%) and less in the other seasons. As summary, it can be indicated that SMP1 and SMP6 are typical winter patterns, while SMP2 is the most frequent summer pattern, with decreasing presence of summer cases and increasing cases of spring cases in SMP5, 4 and 3 respectively. This distribution of cases will be employed later in the interpretation of the results obtained.

3.3. Effect of Saharan Dust Intrusions on MLHs

Regarding the occurrence Saharan dust intrusions, higher MLH is usually observed on those days. This effect can be produced by both the seasonal trend of the episodes and the effect of dust-rich layers over the site. This mixed contribution is difficult to disentangle. As it can be seen in Figure 6, when all cases for the complete four years database are analysed separating by type of day (Clean day or Saharan dust intrusion days) using the methodology explained in section 2.4, a median MLH of 1423 ± 773 m were determined during clean days, and the value increases to 1830 ± 910 m during Saharan dust intrusions (Figure 7a). This increase can be observed systematically in the ABLH, but the behavior of the MLH shows a counter-intuitive feature when the Saharan dust episodes are separated by high and low dust loads. As it can be seen in Figure 7b, high dust load cases, those with ground-level dust concentration higher than 30 µg/m³, as it was explained in section 2.4, produce a lower MLH (1526 ± 876 m) than low dust load cases (1858 ± 914 m). This feature is difficult to explain, so a detailed analysis have been performed.

Firstly, the cases has been analysed only by the identification of Saharan dust days or clean days. Figure 8a shows that Saharan dust days occur in all season, with a pronounced effect over the MLH in autumn, with a mean difference of 427 m. These differences are smaller for winter (315 m), summer (114 m) and spring (6 m). The same study can be perfomed using the SMPs for classification, as shown in Figure 8b. In this case, the main differences are related with the prevalence of some SMPs for occurring in autumn (SMP4: 685 m & SMP3: 558 m), so no advantage is obtained in this season by using the synoptic situations. On the other hand, the differences are very small for one of the winter related cases (SMP1: 20 m) but large for the other (SMP6: 238 m), indicating that the SMP6 pattern is largely affected by Saharan dust intrusions. For the SMP2 cases, the difference is 315 m, so this pattern is also influenced by Saharan dust intrusions.

Figure 9 shows the distribution of clean, low and high dust cases for the six synoptic patterns. As it can be observed, the maximum number of days corresponds to clean days (296 days) for the SMP2, with 108 days of low dust load and 17 days of high dust load for that SMP. For SMP1, no high dust load days occurred in the four years studied, but it happens in SMP6 (12 days), the other prevalent winter case. For SMP3, SMP4 and SMP5, high dust load days were 5, 10 and 6 respectively, therefore, the results obtained using this strategy may be affected by the low High dust cases and the analysis will be more robust if the study is extended to longer period of time.

Taking these into account, the MLH have been analysed for these distribution, obtaining the results shown in Figure 10. In Figure 10a, the distribution by season is plotted, showing that the high dust load days produce lower MLH than both the low dust load and clean days in spring and summer. For the winter cases, the tendency is still increasing with the dust load, although the MLH is the lowest in these cases for the three situations. When the distribution is considered regarding the SMPs, shown in Figure 10b, a pronounced effect is observed for SMP5, occurring mainly in spring, with lower MLH for high dust load cases than both low dust load and clean days. The same effect but less intense is observed for SMP3 and SMP4, and very noticeable for SMP2, occurring mainly in summer, and SMP6, typical of winter and autumn, indicating that dust-rich airmasses behave differently in these synoptic patterns. A possible explanation of this feature is that the aloft dust-rich layers can produce a compression of the MLH, as it has been published by Salvador [45], due to mechanical compression or reduction in radiation reaching the surface.

4. Discussion

In this study, the effect of Saharan dust intrusions on the temporal evolution of the mixing layer heights (MLHs) has been analyzed. Ceilometer profiles were employed to analyze the vertical atmosphere over the site by means of the STRATfinder algorithm to obtain MLH. The dataset was differentiated by different meteorological conditions (seasons, synoptic patterns). Cluster analysis of sea level pressure fields allow the classification of six synoptic meteorological patterns as the main responsible of the weather conditions in the Iberian Peninsula. The mixing layer heights (MLHs) obtained by the STRATfinder algorithm from ceilometer profiles shows clear seasonal effects along the four years analysed (2020 – 2023). A counter-intuitive behavior was obtained, with a systematic reduction of the MLH during the occurrence of strong Saharan dust intrusions, with dust load PM10 concentrations larger than 30 µg/m³. This behavior occurred for all seasons and synoptic meteorological patterns. In general, the distribution of cases by synoptic meteorological patterns provides more insight about the situation, as it allows to differentiate the meteorological conditions suffered by the dust-rich airmass affecting the site atmospheric structure. This effect illustrates the strong direct radiative forcing by desert dust plumes on atmospheric stability within the lower troposphere in particular synoptic situations due to the presence of dust-rich aloft layers. Similar features have been obtained by Salvador [45]. In that study, the impact on air pollutant concentrations and health in Madrid during the period of 2011– 2014 due to variations in the MLH at midday was established. A well correlated decrease between the MLH and a linear increase in the daily number of exceedances of the European Union’s NO2 hourly limit value (200 µg/m³) at hotspot urban traffic monitoring stations was documented. Also an increase in all-natural-cause daily mortality shown a statistically significant relationship with reductions in the MLH. These results are relevant for health advice during Saharan dust intrusion days. It can be concluded that, for a complete description of the temporal evolution of the ML, strongly conditioned by the meteorological conditions, detailed analysis considering synoptic meteorological patterns will yield more clear results. The recent implementation of ceilometer networks will allow for a better characterization of the complexity of ML dynamics at larger scales, offering great potential as a correction tool for MLHs derived from models. This study has also highlighted the need for the use of datasets longer than four years, due to scarcity of cases in several synoptic situations. Future research should focus on characterizing larger geographical areas, taking advantage of the network capacities, in order to improve air pollution dispersion assessments and health alerts to population.