Inventory of African desert dust events in the North-central Iberian Peninsula in 2003-2014 based on Sun photometer-AERONET and particulate mass-EMEP data

A reliable identification of Desert Dust (DD) episodes over North-central Spain is carried out based on AErosol RObotic NETwork (AERONET) columnar aerosol sun-photometer (aerosol optical depth, AOD, and Ångström exponent, α) and European Monitoring and Evaluation Programme (EMEP) surface particulate mass concentration (PMx, x=10, 2.5, and 2.5-10 μm) as main core data. The impact of DD on background aerosol conditions is detectable by means of aerosol load thresholds and complementary information provided by HYSPLIT (Hybrid Single Particle Lagrangian Integrated Trajectory 15 Model) air mass back-trajectories, MODIS (Moderate Resolution Imaging Spectroradiometer) images, forecasting aerosol models, and synoptic maps, which had been carefully reviewed by a human observer for each day included in the DD inventory. This identification method allows the detection of low and moderate DD intrusions and also mixtures of mineral dust with other aerosol types by means of the analysis of α. During the period studied (2003-2014), a total of 152 DD episodes composed of 419 days are identified. Overall, this means ~13 episodes and ~35 days per year with DD intrusion, 20 representing 9.6% days/year. During the identified DD outbreaks, 19 daily exceedances over 50 μg m are reported at the surface. The occurrence of DD event days along the year peaks in March and June with a marked minimum in April and lowest occurrence in winter. A large inter-annual variability is observed showing a statistically significant temporal decreasing trend of ~3 days/year. As a key point, the DD impact on the aerosol climatology is addressed by evaluating the DD contribution to AOD, PM10, PM2.5, and PM2.5-10 obtaining mean values of 0.015 (11.5%), 1.3 μg m (11.8%), 0.55 μg m 25 3 (8.5%) and 0.79 μg m (16.1%), respectively. Almost similar annual cycles of DD contribution are obtained for AOD and PM10 with two maxima, one in summer (0.03 and 2.4 μg m for AOD in June and PM10 in August, respectively) and another in March (0.02 for AOD and 2.2 μg m for PM10), discrepancies occurring only in July and September. It is worth mentioning that the seasonal cycle of DD contribution to AOD does not follow the pattern of the total AOD (near bell shape), meanwhile both PM10 cycles (total and DD contribution) present more similar shapes between them, although a main 30 discrepancy is observed in September. The inter-annual evolution of the DD contribution to AOD and PM10 has evidenced a progressive decrease. This decline in the levels of natural mineral dust aerosols can explain up to the 30% of the total aerosol

Therefore, the contribution of mineral dust in South-Europe is important because of the link between PM 10 exceedances and DD outbreaks.
Once the identification of DD African aerosols is carried out, the next task is to quantify their contribution to the total aerosol load. The evaluation of the contribution of DD episodes to PMx data is viable by means of a chemical speciation analysis (Rodríguez et al., 2001(Rodríguez et al., , 2002(Rodríguez et al., , 2015 but this method requires high man power and has poor temporal sampling . 5 Hence, in order to avoid this expensive technique other methods have been developed using PMx (Escudero et al., 2007;Ganor et al., 2009) and AOD data (Toledano et al., 2007b). As reported by Viana et al., (2010) and taking into account more recent publications (e.g., MAGRAMA, 2013MAGRAMA, , 2015 no more than three methods are currently used with PMx including receptor models (Pey et al., 2013b;Belis et al., 2013). However, these techniques could need to be updated to the measuring site for DD contribution estimates. In a similar way columnar aerosol algorithms can facilitate the apportioning of the 10 different aerosol types (Dubovik et al., 2002;O'Neill et al., 2003) to the total aerosol load.
The advantage of remote sensing techniques, such as sun-photometry, for DD detection is the spectral information recorded by their AOD measurements and given by the Ångström exponent, α. This is a powerful tool in the identification and classification of the different aerosol types (Eck et al., 1999;Toledano et al., 2007a) but also allows "near-real-time" processing of data by means of reasonably sophisticated algorithms (Dubovik et al., 2002;O'Neill et al., 2003) that retrieve 15 aerosol properties. The evolutions of surface PMx and columnar AOD differ in the seasonal cycles (see, e.g., Bennouna et al., 2014;Mateos et al., 2015) and hence their DD impact can also present some discrepancies in that cycle. The PMx sampling used here is based on daily filter records (see, Aas et al., 2013) while sun photometers provide instantaneous measurements of the columnar load but their sampling is limited to daytime cloud-free conditions (Toledano et al., 2007a,b).
The usefulness of a DD inventory is that it opens the possibility of the evaluation of desert dust contribution to the total 20 aerosol load. However, very few studies have accomplished this task over a long period and from a multi-year perspective.
To our knowledge, only the inventory of Toledano et al. (2007b) addressed the DD contribution to AOD between 2000 and 2005 in a Spanish South-western site (El Arenosillo). For PMx data we have found various studies, as those of more recent publication by Salvador et al. (2013; and Pey et al. (2013a). Salvador et al. (2013) reported a DD inventory and the corresponding contribution of DD is determined over the Madrid area over the period 2001-2011, which is extended to 25 several stations covering the whole IP by Salvador et al. (2014). Pey et al. (2013a), with the same methodology as Salvador et al. (2013), analysed the period 2001-2011 for PM 10 at different sites over the whole Mediterranean Basin. In these three mentioned studies, the method used for the evaluation of the DD contribution to PMx is different to that for AOD. Therefore, the development of methodologies for the evaluation of DD contribution is an open area of research, for instance, to obtain the near real-time DD contribution value. 30 Within this framework, the main purpose of this study is to establish an inventory of DD episodes together with the evaluation of their contribution to the total aerosol load, using both AOD and PMx data. This is the first time, to our knowledge, that DD events are identified by the simultaneous use of both columnar (AOD/α) and surface (PMx) aerosol observations. The methodology is applied over one of the cleanest atmospheric areas in south-western Europe, the north-central area ('Castilla y León' region) of the Iberian Peninsula, for a long time period spanning more than one decade (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). Regarding the columnar aerosol data, the reliable measurements performed within the Aerosol Robotic Network, AERONET (Holben et al., 1998), are used. For the study area, the only available aerosol long-term data set is recorded in the Palencia site (41.9º N, 4.5º W, and 750 m a.s.l.). Regarding surface aerosols, we use the high quality particulate matter data recorded by the European Monitoring and Evaluation Programme (EMEP) network in the nearby Peñausende station 5 (41.28ºN, 5.87ºW, and 985 m a.s.l.), with similar background conditions to Palencia. In this way, the use of these two worldwide extended and high-quality networks will ensure the feasibility of implementing the proposed method in other regions. The long-term inventory described hereafter has been employed to establish the main characteristics of the DD episodes in north-central Spain: their climatology, inter-annual behaviour, trends for the number of episodes and associated days, and occurrence under different synoptic scenarios. In addition, the evaluation of the DD contribution to the total mean 10 values of AOD and PMx is also addressed over the period investigated from a climatological and inter-annual perspective, which emphasizes the correlations between both quantities.
Section 2 describes the region of study and the datasets used. The methodology followed in the DD identification and in the evaluation of its contribution to aerosol load is presented in section 3. In subsections 4.1-4.3, the seasonal cycles and interannual evolution of DD events, dusty days, and DD contribution to AOD and PMx are investigated. Subsection 4.4 provides 15 an estimation of the uncertainty of this method and subsection 4.5 describes the synoptic scenarios associated with the arrival of DD episodes to the north-central Iberian Peninsula. Finally, Section 5 summarizes the main findings obtained in this study.
The Palencia site is placed in the autonomous region of "Castilla y León" in the north-central Iberian Peninsula, which is also known as the "Castilian plateau" with an average altitude of ~800 m. This region is the third less-populated community in Spain due to its large area (94,193 km 2 ) and its low population (2.543.413 inhabitants registered in the census in 2012), with a population density just a bit higher than 27 inhabitants per km 2 . Palencia is a small city (100.000 inhabitants) located in the north of "Castilla y León" but the measuring site is located outskirts being surrounded by rural areas, removed from 5 big urban and industrial centres. Hence, this area exhibits an exceptionally clean atmosphere and aerosol observations are representative of the background conditions for the whole region. Therefore, desert dust intrusions can be observed since they significantly modify the background aerosol properties.
In order to fill the gaps in the Palencia AOD database, the site "Autilla" (42.00ºN, 4.60ºW, and 873 m a.s.l.) close to Palencia (7 km apart) has been used (for details see Bennouna et al., 2013). This site is used by GOA (Grupo de Optica 10 Atmosférica) as the calibration platform for Cimel sun photometers within the AERONET-EUROPE infrastructure and also works as a routine measurement site. Under these considerations, the columnar aerosol data series used in this study is consistent and allows one to perform the inventory of DD events in this region.

EMEP network and PM database 15
Daily PM 10 and PM 2.5 measurements provided by the EMEP network constitute the second core database used to carry out this study (see Figure 1 and Table 1). This network has the objective of regularly providing qualified scientific data to interested organizations in order to analyse and assess the transboundary transport and emission of pollutants (e.g., Aas et al., with their atmospheric and background conditions make possible the joint discrimination and evaluation of these observations of AOD and PMx for the detection of DD intrusions. Furthermore, an analysis about the air masses at Peñausende and Palencia sites (not shown here) has been carried out, corroborating that the geographical distance between them is negligible for the analysis of regional quantities such as AOD, water vapour, and ozone column, among others.
The use of AOD 440nm , α, PM 10 , and PM 2.5 observations provides a comprehensive database to carry out an analysis of aerosol 5 load and particle size, both at the surface and in the whole atmospheric column. Table 1 presents a detailed description of the number of days with available data every year of each database. Overall, PMx presents a larger amount of days with data than AOD. Particularly, the year 2003 presents the lowest AOD sampling (42% of 365 days) because of certain gaps just at the start of the sun photometer measurements.

Ancillary information 10
To carry out a more accurate evaluation and discrimination of days that constitute a desert dust intrusion, ancillary information is also considered. Air mass back-trajectories arriving at Palencia at 12.00 UTC have been calculated with the HYSPLIT model (Hybrid Single-Particle Lagrangian Integrated Trajectory), version 4 (Draxler et al., 2014;Stein et al., 2015). Due to the fact that desert dust aerosols can be transported to altitude levels higher than the boundary layer, backtrajectories have been calculated for three heights (500, 1500 and 3000 m above ground level) and analysed 5 days back in 15 time (120 h), using the model vertical velocity in the calculations. The meteorological database used as input for HYSPLIT is the Global Data Assimilation System (GDAS) GDAS1 dataset (e.g., Su et al., 2015). These three levels are commonly used in these studies to represent the air masses near the surface, in the boundary layer and in the free troposphere, in order to follow the vertical transport of aerosols.
Information about cloudiness is obtained from MODIS (Moderate Resolution Imaging Radiometer) rapid response imagery 20 products (https://earthdata.nasa.gov/data/near-real-time-data/rapid-response). In addition, GIOVANNI (Geospatial   Interactive  Online  Visualization  ANd  aNalysis  Infrastructure) MODIS AOD aerosol maps (http://giovanni.gsfc.nasa.gov/giovanni/) and those provided by the AERONET website (http://aeronet.gsfc.nasa.gov/cgibin/bamgomas_interactive) are used to determine the extension and path followed by the mineral dust air masses for a DD identified event. The NAAPS Global Aerosol model (Navy Aerosol Analysis and Prediction System; available at 25 http://www.nrlmry.navy.mil/aerosol/) is also used to corroborate if model forecasts also detect a given DD episode over the study area.
As shown in previous studies (Escudero et al. 2005;Toledano et al. 2007b) desert dust intrusions over Spain take place under certain synoptic scenarios (see Section 4.5 for further details). Through the Earth System Research Laboratory of NOAA (National Oceanic and Atmospheric Administration), the plots of the geopotential height at 700 hPa and mean sea level 30 pressure are obtain in order to evaluate the synoptic scenario associated for a given DD episode among the four possibilities (see Section 4.5) proposed by Escudero et al. (2005).

Detection of desert dust episodes
This study is based on instantaneous AOD 440nm and α values, as well as daily PM 10 and PM 10 /PM 2.5 ratio data. The method for the detection of Desert Dust (DD) intrusions is a manual inspection of the evolution of these four quantities together with the origin of the air masses at the three levels of altitude at 500, 1500, and 3000 m a.s.l. and the auxiliary material of AOD 5 MODIS maps, aerosol models and synoptic scenarios. The methodology for detection is similar to that applied in Toledano et al., (2007b) with the added information of PMx data, and not so much different from that used in other studies also based on a set of different observations (Escudero et al., 2005;Pace et al., 2006;Kalivitis et al., 2007;MAGRAMA, 2013). The difference between these methods lies on the weight played by each quantity, and the way to analyse the information. For example in our case the AOD-PMx data is the first and primary information, but in other studies the key variable is the 10 origin of the air masses (Pace et al., 2006;Valenzuela et al., 2012). Meteorological products and forecast aerosol models can also be used for this task (MAGRAMA, 2013). Although automatic methods can be applied in the DD identification, a visual inspection should be performed to corroborate each classification.
This study has been carried out as a year by year service to the "Consejería de Medio Ambiente" of the Autonomous Community of 'Castilla y León' by means of two Research Programmes from 2006-2013 about the "Discrimination, 15 characterization and evaluation of desert dust outbreaks over 'Castilla y León' region". These programmes aim to help implement the "Environmental Quality Improvement Policy" of the EU by the National and Regional Governments of Spain.
The experience gained with this year by year identification provides the final DD inventory presented in this study.
Certain thresholds have been established to identify those conditions which stand above the clean background over the study area. Hence, the choice of these thresholds is based on the aerosol climatology of the site from our investigations (Bennouna 20 et al., 2014;Mateos et al., 2015) and previous results (see, e.g., Querol et al., 2009;. The mean values for the longterm period 2003-2014 are 0.13 ± 0.09 for AOD 440nm and 10 ± 9 µg m -3 for PM 10 . Hence, to detect the DD intrusions a visual inspection of the entire database is performed. When a day shows a group of number of points of the instantaneous AOD ≥ 0.18 and/or the daily PM 10 ≥ 13 µg m -3 , that day is further investigated. The AOD threshold corresponds to the mean value plus half of the standard deviation, approximately, meanwhile that for PM 10 is the mean plus one third of the standard 25 deviation. Hence, these thresholds must be taken as "warning flags" in the sense that these values alone do not define the classification as a dusty day. They also need the ancillary information given by the air mass backwards trajectories, satellite images, weather maps, and model forecasts to determine and corroborate the origin of aerosols and synoptic conditions. Therefore, with all this information the human observer decides if a day must be included or not in the DD inventory. In parallel to the above analysis, the evolution of the α quantity is also checked, allowing the identification of two different 30 types of DD intrusions. Those days displaying α ≤ 1.0 in most of the instantaneous columnar data are identified as the "purer" desert dust intrusions and they are denoted by a D flag. Those days with α values in the interval 1.0 ≤ α ≤ 1.5 (which have been classified as dusty days by their aerosol load and/or the ancillary information) present a mixture with other types (e.g., clean continental aerosols) and they are denoted by a MD (Mixed Desert) flag. The MD event days may either be a part of an intense event (generally of D type) or form by themselves a low-moderate intensity event. Previous studies have also shown that DD intrusions in the Mediterranean Basin can present moderate AOD associated with large α values (e.g., Pace et al., 2006;Tafuro et al., 2006;Boselli et al., 2012). The limit of α ≤ 1 for identifying coarse particles has been established by previous studies in different worldwide areas (e.g., Eck et al., 1999;Dubovik et al., 2002, Meloni et al., 2007Boselli 5 et al., 2012), making this threshold suitable for our study area. It is important to emphasize here that the α parameter allows a more accurate identification of desert dust events, mainly those of low intensity, having less of the characteristics identifying a desert origin (overall, the larger the AOD of a desert event, the lower the α), generally mixed with other aerosol types which are not accounted for in many DD studies (e.g., Gkikas et al., 2013). These DD events of low intensity are more difficult to detect because of the low signal of mineral dust aerosols and hence it is more difficult to evaluate their 10 contribution or impact on AOD and PMx daily values. In fact, although the aerosol load threshold used in this study might seem very low, we must note again here that all these quantities are manually supervised by the expert-observer who will take the final decision of the inclusion or not of a dusty day, bearing in mind all the available information given by these data and all the complementary material.

15
Another problem to be considered in the identification of a dusty day is when AOD and PM 10 measurements show different information. For instance, AOD fingerprints seem to indicate possible desert dust presence meanwhile PM 10 not. It must be taken into account that PM 10 quantity does not necessarily follow the same temporal behaviour as the AOD and possible time delays in PM 10 concentration due to deposition processes can occur. Desert dust events can reach the IP at high altitude layers (e.g., above 2000 m, Gkikas et al., 2015), and dry deposition can last various hours/days. Assuming an average speed 20 of deposition of around 0.6 cm s-1 (Zender et al., 2003), DD particles can remain in the troposphere up to two days after an episode ends (Escudero et al., 2007). Hence, these possible day delays between columnar air masses and surface PM 10 are very variable (e.g., Kalivitis et al., 2007;Pey et al., 2013a). Therefore, AOD and PM 10 observations must be considered as complementary information to detect mineral dust aerosols since with the evolution of both quantities it is easier to identify DD fingerprints. 25 It is worth mentioning that in the general detection of both intensity and duration of a particular event, DD causes an increase in the AOD 440nm and PM 10 values, which then surpass or just reach the corresponding threshold, together with a decrease in the α and PM 2.5 /PM 10 values, giving rise to an increase of the mean size of the particle distribution. The duration of each intrusion can be established, since the low typical values, characteristics of the regional background, are recovered when the event finishes. The central or most intense days of each event are easy to detect due to the large increment of aerosol load 30 with large particles (α even close to zero) but low to moderate events are more difficult to detect. Although the events vary in their nature, the first and last days of a DD event show a low or moderate signature of mineral dust particles because of its mixture with other aerosol types (e.g., clean continental aerosol in our study area), with the exception of the strong DD events which generally have a notable impact on the aerosol load levels from the first day of the episode.
In the inspection of the instantaneous columnar dataset, non-reliable records are identified and removed due to their high spread in the data likely attributed to cloudy conditions. In the corroboration of some critical decisions made by the expertobserver the ancillary information (see section 2.3) constitutes a key point of the methodology. For instance, the verification of cloudy conditions can be supported by means of a comparison between AOD instantaneous data from different levels, 1.0 (all-sky conditions) and 1.5 (cloud-screened data), and the visualization of cloud systems in MODIS true colour and cloud 5 product images. If there are signs of cloud presence, instantaneous AOD data are carefully checked to discern between cloud-contaminated data and a DD intrusion.
Once aerosol load measurements for a given day indicate the likely classification as a dusty day, the air mass backtrajectories (calculated as described in Section 2.3) are visualized to check if the origin of the path or the path followed crosses the north-African region and/or its surroundings. Therefore, the air mass back-trajectory analysis and the 10 geopotential maps (establishing a particular synoptic scenario, see Sections 2.3 and 4.5) lead to the final decision with respect to a DD event day classification, even in those days showing cloudiness. Finally, to help in the understanding of the general situation about the geographical distribution of the aerosol plumes, AOD MODIS maps and NAAPS forecasts are also visualized for ensuring the final choice. The consistency of the information used provides a reliable identification of the DD event. 15 It is worth mentioning here that the final decision to include a day as D or MD is made by the human-observer with all the available information at hand. Perhaps this methodology is not the most adequate to apply for a big area with a high number of sites and long-term databases, but it is necessary for developing methodologies, because it will allow validation of other more automatic methods (e.g., those only based on threshold criteria), most of them using satellite observations (e.g., Gkikas et al., 2013Gkikas et al., , 2015. 20

Evaluation of desert dust contribution to total AOD and total PM 10 concentration
Once the DD inventory is established, the evaluation of the DD contribution can be addressed at seasonal and annual scales.
Following Toledano et al., (2007b), the contribution of the DD events to AOD can be obtained as the difference of the multiannual monthly means considering all days and the corresponding value without including the desert dust cases. This 25 procedure was also used for PM 10 data in a 3-year evaluation of net DD contribution in several sites in the Mediterranean Basin (Querol et al., 2009). In this study, the annual cycle of the DD contribution to AOD/PMx is evaluated with this same methodology over the entire period 2003-2014, using the DD event days classified in the inventory. Furthermore, the relative DD contribution to AOD/PMx can be obtained by dividing each one with respect to the corresponding total AOD/PMx value. Regarding the seasonal evaluations, the classification is as follows: winter (Dec-Jan-Feb), spring (Mar-Apr-May), 30 summer (Jun-Jul-Aug), and autumn (Sep-Oct-Nov). Analogously, the yearly AOD/PMx means excluding dusty days are subtracted from the yearly AOD/PMx means for all days to obtain the DD contribution at annual time scale. This method assumes the entire daily aerosol load (both surface and columnar) due to DD aerosols, being the contribution of regional background aerosols also included. Thus, suitable time scales for this kind of DD contribution calculation are the annual and multi-annual monthly means over 12 years of data. Those evaluations for every single day or month can be addressed using other methods, for instance the determination of percentile 40 of the time series without dusty days to evaluate the background conditions that are subtracted from PMx levels (Escudero et al., 2007). This method has been taken 5 as the standard by the European Commission for the evaluation of DD contribution to aerosol load at the surface (Viana et al., 2014;MAGRAMA, 2013MAGRAMA, , 2015. The methodology used in this study leads to lower uncertainty in the annual cycle evaluation since there is good data coverage for the multi-annual monthly sampling (12 years). However, in the evaluation of year-by-year DD contribution to aerosol load, a given year can present a low coverage leading to a higher uncertainty. The inclusion or not of an uncertain 10 DD event (e.g., data contaminated with clouds where cloud optical depth is assigned to aerosol AOD) can substantially modify the corresponding yearly mean as it has been shown in Bennouna et al. (2014). This source of uncertainty must be considered in the temporal trend evaluation. It is not easy to establish an adequate methodology to evaluate the DD aerosol contribution to aerosol load (Viana et al., 2010) and much less with its corresponding associated error. A further investigation is necessary about this subject. A discussion about the uncertainty of our approaches in the DD identification 15 and in the evaluation of DD contribution to aerosol load can be found in subsection 4.4.

Evaluation of the number of episodes and dusty days
The inventory of desert dust intrusions includes: information on each episode and its associated days; the daily mean AOD, 20 α, PM 10 , and PM 2.5 ; cloudiness, synoptic scenarios, and air mass origin at three altitude levels (500, 1500, and 3000 m a.s.l.).
Tables 1 and 2 show the information used to classify DD events and the main statistics for this inventory, respectively.
The PM 10 sampling presents the best coverage of the measuring time period with 93.1% of the days, AOD is available 67.2% of the time, and the coincident sampling is available 63.2% of the time. As can be deduced from Table 1, the majority (51.2%) of the DD event days composing the inventory are noticeable in both AOD and PM10 datasets. However, 46.4% of 25 the total detected days are over the required thresholds only in one quantity (AOD or PM10). This is the advantage of the proposed inventory because DD outbreaks are identified with two complementary quantities about the aerosol load. The reasons behind this 46.4% (19.1% only with AOD and 27.3% only with PM10) are due to time delays between columnar and surface levels related to deposition phenomena and also are due to the lack of AOD or PM10 measurements. Therefore, if the inventory is addressed by only one quantity, a large number of dusty days can be lost in their identification. Finally, a 12 smaller number of cases (2.4%) are identified as dusty days using the ancillary information when AOD and PM10 data are not available.
The smaller amount of available data in the AOD time series (see Table 1), in comparison with PM10, is not a major handicap in DD detection. There are several years (2004, 2007, 2008, 2010, 2011, and 2012) with more DD days detected only by AOD than only by PM 10 , in spite of the smaller AOD sampling (between 53 and 103 days less per year). However, 5 years 2003 and 2006, which present less than 200 daily AOD data, require the use of PM 10 to better identify DD intrusions.
As reported in Table 2, during 2003-2014, a total number of 152 episodes have been identified, composed of 418 days.
Among them, 242 days have been classified as days with desert aerosols (D) and 176 days with mixed aerosols (MD).
Overall, this means 13 episodes and 35 days per year with desert dust intrusion, representing 9.5% of the days each year. The duration of DD episodes is very variable, ranging from 1 to 13 days, but a value of 2.7 days is obtained as the mean episode 10 duration. Due to the high variability in the intensity of these intrusions it is difficult to distinguish when an event has ended or its intensity has simply fallen below our threshold. During summer, the recirculation of air masses in the IP is very frequent and the DD episodes are subject to large intensity variations. We have considered separate DD episodes when there is, at least, one day that does not meet the DD requirements between two DD episodes.
Our percentage of dusty days of 9.5% is lower than that reported by Salvador et al. (2013), which is around 18% ( Figure 2b is similar in shape to that reported by Salvador et al., (2013) with the exception of September and also similar to that of Escudero et al. (2005) with the exception of October. Concerning the two types of DD conditions 30 distinguished in our inventory, D type controls the annual cycle in the March maximum and April minimum, while MD controls the evolution between August and October. Some features mentioned above regarding the seasonal behaviour of DD events for the north-central Spanish region are also observed for other areas of the IP. For instance, the March maximum and April minimum are common features in southwestern (Toledano et al., 2007b;Obregón et al., 2012), north-eastern (Escudero et al., 2005Papayannis et al., 2008), and central (Salvador et al., 2013) Spain. In spite of the different measuring instruments (sun-photometer, lidar, or particulate mass concentration), time periods, and methodologies employed in the DD event identification, its impact over (almost) all 5 the Iberian Peninsula follows the same pattern with the two maxima, one in late-winter/early-spring (March) and the other in summer, and the accentuated minimum of winter. Several minor discrepancies are found for the rest of the year; for example, the maximum number of DD event days during summer months with a local minimum (in July) between them seems to be a characteristic of the annual cycles of the south-western (Toledano et al., 2007b) and north-central areas. Conversely, the eastern region presents a local maximum in October (Escudero et al., 2005;Pey et al., 2013a) which is not seen in our study 10 region. These results confirm that different areas have different aerosol properties in the IP (Mateos et al., 2015). The larger occurrence of dusty days in summer and in certain spring and fall months is also observed in other Mediterranean areas using lidar networks to identify the occurrence of DD outbreaks (e.g., Mona et al., 2006;Papayannis et al., 2008).

Inter-annual variability and trends of the number of episodes and days with DD
The year-to-year variability for both number of episodes and days is reported in Table 2  To quantify the decreasing trends in the number of episodes and associated days, the Theil-Sen estimator and Mann-Kendall test for significance have been used. The trends for the number of DD episodes and days are reported in Table 3 for the 30 yearly values. A statistically significant trend at the 95% significance level presents a p-value below 0.05 (e.g., Sanchez-episodes per year with a p-value around 0.03 (~97% of significance level). These figures corroborate a significant decrease in the DD events seen in the north-central area of the Iberian Peninsula over the past decade. This result is in line with the findings obtained by Gkikas et al. (2013) for the whole Mediterranean Basin using MODIS data between 2000 and 2007 and considering only very intense DD events selected by a high AOD threshold. Figure 4 and Table 4 show the annual cycle of the DD contribution (small red bars in the figure) together with the multiannual monthly means considering all days and only days without DD aerosols (the difference between these two values gives the DD contribution). Overall, the mean DD contribution to AOD is 0.015 or 11.5% in 2003-2014.

Annual seasonal cycle
The total AOD annual cycle representing the climatology follows the well-known pattern previously reported and explained 10 for the Palencia site (see, e.g., Bennouna et al., 2013;Mateos et al., 2014a). To summarize: the increasing values from January to June (just when the maximum is found) with a slight reduction in May and a decreasing trend to the end of the year, provide almost a well-defined bell shape. Concerning the climatology with the DD episodes excluded, it preserves the bell pattern found before for the general case, except for some minor discrepancies. For instance, the change between May and June is not noticeable for the curve with the DD excluded, in contrast with the larger increment observed for the general 15 case.
However, the seasonal pattern followed by the DD contribution to AOD is considerably different to these two latter quantities. Two maxima are observed during the annual cycle: the first one in March (late winter/early spring) with 0.018 or 13.4% and the strongest one occurring in summer period (June and August), ~0.027 or ~17%. Together with these maxima, there are two local minima: in April-May (around 0.014 or 9.5%) and in July (0.018 or 12%). After August, a progressive 20 decline of the DD contribution is observed with the minimum in winter (December and January show similar values about 0.004 or 5.4%). It is worth mentioning here the different characters of the two local minima occurring in April-May and July, the former is more general of the IP (linked to the precipitation cycle) while the latter is more typical of the central and south-western areas of Spain. For instance, the July minimum seems to be related with the arrival of drier air masses in the low troposphere as it is observed in the precipitable water vapour cycle (Ortiz de Galisteo et al., 2013). 25 The annual cycle of the DD contribution to the Palencia site (representing north-central Spain) presents a similar shape to that obtained in the "El Arenosillo" site (south-western area) by Toledano et al., (2007b) for an inventory of 6 years, from 2000 to 2005. This is an important result in two aspects, one related with the shape of the annual cycle or seasonal behaviour and the other one related with the different contribution of north and south areas of the IP. In relation to the geographical gradient of African dust contribution, a quantitative difference is observed between these two areas. The total AOD signal is 30 clearly impacted by DD events in the southern Iberian coast (with relative contributions being over the 30%), while in the north-central region the DD influence is weaker, thus a south-north decreasing gradient over the IP is observed regarding the DD contribution to AOD values. This behaviour is well known in the IP by earlier aerosol studies based on PMx data (Querol et al., 2009;Pey et al., 2013a;Salvador et al., 2013Salvador et al., , 2014 but this is the first time this is confirmed by an inventory of AOD data.

Inter-annual variability and trends
With respect to the inter-annual change of the DD contribution to AOD, Figure 5 and Table 5  The temporal trends in total AOD and in the DD contribution to AOD are also evaluated and shown in Table 3. The decrease of the total AOD in the Palencia site in 2003-2014 is 0.006 AOD-units per year (with a p-value <0.01) or 4.6% per year, which is in line with previous findings for the same site by Bennouna et al., (2014) and Mateos et al., (2014b) for shorter (4 and 3 years, respectively) periods. With respect to the DD contribution to AOD, a decrease of 0.0019 AOD-unit per year (p-25 value = 0.02) or 11.2% per year is calculated. Therefore, this rate represents the 30% of the total AOD decreasing trend.
Hence, the natural decrease of DD aerosols has notably affected AOD levels over north-central Iberian Peninsula during the study period.

Annual seasonal cycle
In the same way as for AOD, the contribution of desert dust events to mean values of PM 10 , PM 2.5 , and PM 2.5-10 have also been calculated. The annual cycle and the inter-annual evolution of these three quantities and the corresponding DD contributions are reported in Tables 4 and 5 and also illustrated by Figures 6 and 7, respectively. 5 The DD contribution to the total PM 10 , PM 2.5 , and PM 2.5-10 is not usually evaluated at the same time. To our knowledge, this is the first time that fine and coarse mode contributions are evaluated in a long-term desert dust inventory of this type. Furthermore, the temporal trends for the inter-annual DD contributions are also discussed. It is worth mentioning here that as PM 10 and PM 2.5 are obtained from independent filters (see Section 2.2) while PM 2.5-10 is only available with simultaneous PMx data, the amount of data used in the evaluation of DD contribution for each quantity slightly differs. 10 According to Table 4, the mean DD contributions to PMx during the study period are 1.3 µg m -3 (12%) for PM 10  April. This general behaviour for the entire dataset is also followed if the DD events are excluded. The evolution of these two latter curves is also followed by the DD contribution to PM 10 . The largest DD contribution is observed in March (2.2 µg m -3 or 20%) and summer months, June to August (~2.3 µg m -3 or ~17%). The months of April and May (~0.9 µg m -3 or ~9%) display a notable decrease with respect to March. After summer, there is a sharp fall in September (1.2 µg m -3 or 10%) producing a local minimum, and beyond October a progressive decline leading to the weakest effect (<8%) of the African 25 intrusions during winter months (DJF). The maximum relative DD contribution to PM 10 can reach 20%, which is within the range (10%-50%) observed by Pey et al. (2013a) for the eastern Spanish coast. Comparing the seasonal cycles of DD contribution to PM 10 in the latter area with respect to north-central Iberian Peninsula, some common features appear (March maximum, April/May decrease, summer increase, and September drop) but the maximum in October seen in the Mediterranean coast does not happen in the north-central area. 30 Even though both AOD and PM 10 express the aerosol load, these quantities present noticeable differences. To facilitate the comparison of the results shown above, Figure S1 (supplementary information) shows together the annual cycles of AOD and PM 10 total means and their DD contributions. The annual cycle of the two quantities, total AOD and PM 10 , for the complete dataset follows a similar behaviour between August and March, with differences in April (local PM-minimum) and May (local AOD-minimum) , and a slight different evolution in June-July. These discrepancies between these quantities lead to a moderate-high correlation coefficient of 0.82 between AOD and PM 10 , but their physical meaning is uncertain taking into account the discrepancies in the two annual cycles. The seasonal cycles of DD contribution (both absolute and relative), to AOD and PM 10 differ in July (with a local minimum of AOD) and September (sharp fall of PM 10 ). Furthermore, it is 5 worth mentioning to note that the maximum of March is more intense for the DD contribution to PM 10 than to AOD. Hence, the correlation factors between the DD contribution to PM 10 and to AOD are moderate-high: 0.84 and 0.74 for the absolute and relative values, respectively.
The fine mode, represented by the PM 2.5 data, follows the same pattern as PM 10 in the total and DD contribution curves (Table 4 and Figure 6b). The DD contribution to PM 2.5 is below 10% for most of the year, with a mean value of ~9%. 10 The total coarse mode (PM 2.5-10 ) curve is also similar to that obtained for the total PM 10 , although the mean contribution of the DD events is 16% of the total PM 2.5-10 , faced to the 12% of the PM 10 (see Table 4). The DD contribution to PM 2.5-10 (Table 4 and Figure 6c) exhibits a strong maximum in March (1.7 µg m -3 or 33%), a reduction in April and May (around 14%), large values in June (1.4 µg m -3 or 25%) followed by a weak decrease in July and August (1.3 µg m -3 or 21%), and low values in autumn and winter. 15

Inter-annual variability and trends
The inter-annual variations of total PM 10 , PM 2.5 and PM 2.5-10 and the corresponding DD contributions to these PMx concentrations are plotted in Figure 7 and reported in Table  during that year (see Table 1). So far, no reasonable explanation has been found for the strong fall between 2004 and 2005 in 30 the DD contribution to AOD despite the fact that total AOD and PM 10 display the same behaviour. The DD contribution to PM 10 is notably larger than that obtained for AOD in 2014. The high inter-annual variability of these quantities highlights the necessity of longer time periods to assess this kind of relationships, but bearing in mind that the net contribution of DD aerosols is represented by very low values with a high uncertainty, hence this variability falls within the expected range of change. These results are of interest for long-term studies of columnar and surface aerosol loads in relation to their evolution and trends for climate studies because tropospheric aerosols have a strong regional signature and the area studied presents exceptional background conditions representative of the Western Mediterranean Basin.
The weak impact of the DD events on the PM 2.5 levels (fine mode, see Figure 7b and Table 5) is reflected in the low relative 5 contribution with only three years (2003, 2005, and 2006) presenting values higher than 12%. The last years of the period analysed (2009-2014) present a notable low DD contribution to PM 2.5 below 7%. On the contrary, PM 2.5-10 ( Figure 7c and Table 5) still follows the PM 10 pattern. The initial years are the ones with the largest contributions (around 27% until 2006) while 2013 shows the minimum values (around 5%) together with 2009 (~7%).
There is a decreasing trend of all the quantities shown in Figure 7. The general decrease of PMx levels has been previously 10 reported for the Peñausende site and for shorter periods (e.g., Barmpadimos et al., 2012;Bennouna et al., 2014;Mateos et al., 2014b; and it has been corroborated with the temporal trends obtained in this study (see Table 3). Regarding the DD contribution, the fall in the three quantities is quantified as around -10% per year. In particular, the DD contribution to PM 10 has decreased by an absolute amount of 0.14 µg m -3 per year (p-value of 0.06) and 0.08 µg m -3 per year (p-value < 0.01) for PM 2.5 . The reduction observed in the DD event days (see subsection 4.1.3) has also led to a significant fall of the total particulate matter. Comparing the temporal trends of the PM 10 DD contribution and the rate for the total 20 quantity, the DD impact has caused 30% of the total PM 10 decrease in north-central Spain. This percentage is smaller (about 21%) for the PM 2.5 case. In the north-eastern region,  showed that crustal matter accounted for 14% of the total PM 2.5 decrease between 2001 and 2012.

Estimation of associated uncertainty of the methodology
No quantification has been done about the associated uncertainties in the number of events and associated days in most of 25 the reported bibliography. The same happens for the uncertainty linked to the DD contribution, which can be evaluated as a consequence of the earlier error of DD detection. One way to estimate a possible range of the real uncertainty is the comparison of results obtained from different methodologies. This task was addressed by the above mentioned study of Viana et al. (2010) showing relative differences in the number of dusty days about 12% and 28-50% for the DD absolute contribution. A big step took place when the proposed methodology by Escudero et al., (2007) was taken as the official 30 standard method. However, the 30 days moving percentile used to establish the regional background has been changed from 30% (reported by Escudero et al., 2007) to 40% (Pey et al., 2013a;Salvador et al., 2013Salvador et al., , 2014. There is evidence that this percentile may be site dependent thus demonstrating the difficulty of this evaluation. Otherwise, it must be borne in mind that a big difference exists between the Escudero et al., (2007) methodology and that applied by us. This subsection describes a first attempt to estimate the uncertainty associated with the method used in our study.
Fingerprints of each DD event day are visible on at least one of the quantities related to aerosol load (columnar or surface) analysed in the inventory evaluation (see Section 3), plus the additional informational of air mass back-trajectories, satellite images, and synoptic scenarios. Usually, several of these variables simultaneously corroborate the DD presence, especially 5 due to the low background values that characterize the north-central Spanish region. Therefore, the thorough inspection of all the information provided by different sources at the same time causes the error in the DD identification to be minimal. From our experience during these 12 years of data, we consider that possible error sources can be, mainly, the following: gaps in the data series, classification or not of a day when the aerosol load is close to the threshold values, and uncertainty of the instrumental techniques and the ancillary tools. Therefore, we can estimate that about 3-5 days per year could be missed in 10 the annual sum of dusty days. This assumption is based on our long-term expertise in this evaluation when the DD inventory has been re-evaluated to ensure its accuracy. So the associated relative uncertainty, considering the average of 35 DD event days per year, is ~9-15%, which is in line with the results reported by Viana et al. (2010) as mentioned above. This estimation gives a realistic range for the error associated with this methodology of visual inspection. The 5 days per year uncertainty (or 15%) can overestimate the real error, but even this percentage can be considered as acceptable as the 15 maximum average error. Regarding the sum of dusty days in the seasonal cycle, the same range of error can be assumed in every monthly inter-annual value.
The possibility of missing these few days with DD fingerprints (~3-5 per year and per inter-annual month) leads to an uncertainty in the evaluation of the DD contribution to AOD values. Hence, to quantify the uncertainty in the seasonal cycle of the DD contribution to AOD each inter-annual monthly database is extended adding 9% of DD event days (considering 20 that 3 days in are missed). For these "extra" days the AOD is assumed as the mean value during the DD events in that month.
For instance, four days are added in June with a mean AOD of 0.27 and one day is added in January with AOD 440nm = 0.18.
The DD contribution is calculated for the new data series, evaluating the differences with those values shown in Section 4.2 (from the original database). The results show a small change in the DD contribution to AOD, always below 0.002. For instance, for June the relative uncertainty caused by the added days is 6.7% (the absolute DD contribution for the original 25 evaluation is 0.027). However, those months with less absolute DD contribution to AOD cause a relative difference between 15% and 20% (such as January and December). Overall, the mean uncertainty is 0.0013 or 9.7% when the uncertainties of the twelve multi-annual monthly values are averaged. This relative uncertainty is in line with the 10% calculated by receptor modelling studies (Viana et al., 2010). The same procedure is applied for the inter-annual DD contribution to AOD. On average, the inclusion of 9% DD extra days causes an uncertainty of 0.0014 or 8.3%. If the assumption of missing 3 days per 30 year is even increased to 5 days per year, the uncertainties caused on the DD contribution to AOD values only increase up to 14%. Hence, the reliability of the method followed here is demonstrated because of the low changes of the results when the DD inventory is added with possible dusty days missed.
In the same way, the study of the uncertainties of the DD contribution to PM 10 is also addressed with the same method (adding 9-15% extra DD event days). The results for PM 10 indicate a mean uncertainty of 0.1-0.13 μg m -3 or 8-14% in the evaluation of both annual cycle and inter-annual evolution. Hence, this relative uncertainty can be also extrapolated to the PM 2.5 and PM 2.5-10 DD contributions proving the feasibility of this method.

Analysis of the synoptic scenarios during desert dust episodes 5
Using the ancillary information used in the final choice of the DD identification, the synoptic scenarios that favour the arrival of air masses originating in the north of Africa are also studied. These scenarios are those defined and described by Escudero et al. (2005): via the Atlantic arch (North Africa High Located at Surface Level, NAH-S), directly from north-Africa by a deep low pressure (Atlantic Depression, AD) or by a convective system (North African High located at upper levels, NAH-A), and from the Mediterranean area (North African Depression, NAD). Overall, the geographical positions and 10 heights of the high and low pressure systems produce the mineral aerosols to reach the IP. Figure 8 presents the annual cycle and inter-annual variability of the number of episodes associated with each synoptic scenario. The synoptic scenario of each episode has been established considering all the daily meteorological maps during the episode.
The synoptic scenario analysis of the DD events (see Figure 8a) has shown a predominance of the NAH-A (81 out of 152 episodes), in particular, during the warm season (from May to October). This scenario is produced by an intense solar 15 heating of the Saharan desert. These air masses present large DD loads which can arrive at high altitudes (up to 5 km a.s.l.).
In our study region, the NAH-S scenario governs (38 out of 152 episodes) the DD intrusions between December and April (being also significant in October) and generates transport at the lower atmospheric levels (generally below 1 km a.s.l.). The AD scenario plays a minor role (24 out of 152 episodes) but with an influence confined between February and May, September, and November. The NAD scenario only presents an important contribution in March and December (9 out of 20 152 episodes).
The fingerprints of the evolution of these synoptic scenarios are reflected in the climatology of the DD episodes shown in Figure 2. The rapid increase in DD events in March (see Figure 2) is caused by a larger influence of NAH-S (3 to 5 DD events with respect to February), the marked appearance of NAD (3 events), and a slight increase of AD (2 to 3 DD events with respect to February). The synoptic situation in April changes and the NAD scenario almost disappears while NAH-S 25 and AD increase their influence. The local summer minimum in July is caused by the lower occurrence of the NAH-A conditions. Previous studies have found this minimum for other columnar quantities, such as the vertical precipitable water vapour (Ortiz de Galisteo et al., 2013). The absolute DD event minimum of November is caused by the total disappearance of the NAH-A scenario.
Comparing these results with previous inventories performed in other geographical areas of the IP, the synoptic scenario 30 climatology presents some discrepancies. Toledano et al. (2007b) have also found for the "El Arenosillo" site (south-western IP) in the period 2000-2005 a predominance of the NAH-A conditions during summer. However, the role played by the NAH-S seems to be minor during winter compared to the north-central area. The DD inventory in the north-Mediterranean Spanish coast has been analysed by Escudero et al. (2005Escudero et al. ( ) between 1996Escudero et al. ( and 2002. They also obtained the major predominance of the NAH-A during summer, although the NAD scenario shows a notable impact on the DD events in May and November. These outbreaks arriving from the Mediterranean area are also reported in the months of February, March, and November in the "El Arenosillo" inventory. Inter-annual distribution of DD events and the four synoptic scenarios (see Figure 8b)

Conclusions
In this study, a methodology to obtain a reliable identification of DD intrusions is proposed and applied to the north-central area of the Iberian Peninsula. Long-term datasets of AOD and PMx for the background sites of Palencia and Peñausende (representative of the study area) have been used as core information for the detection of desert dust intrusions in this area during an 12-year period (from January 2003 to December 2014). The analysis of ancillary information, such as air mass 20 back-trajectories at three altitude levels (500, 1500 and 3000 m a.s.l.), MODIS-AOD and true colour images, and meteorological maps, has been used to establish the duration of each desert dust episode, creating a reliable inventory of desert dust episodes. The main conclusions can be summarized as follows: 1. The simultaneous consideration of surface and columnar aerosols has been shown to be a reliable tool in the DD identification. More than a half of the inventory has been detected by AOD 440nm and PM 10 data at the same time. 25 However, each quantity can provide DD detection by itself in a large number of cases (114 and 80 out of 418 days detected by only PM 10 and AOD data, respectively). The smaller coverage of AOD sampling is not a major handicap in this process. days in April, a notable increment between May and September and a progressive decline to the absolute minimum in winter, with the absolute maximum in June and local minimum in July/August. Inter-annual variability of the number of DD episodes and dusty days is high, ranging between 7 episodes (15 dusty days) in 2013 and 17 episodes (67 dusty days) in 2006. A temporal trend of -2.7 dusty days per year (95% significance level) points out its decrease between 2003 and 2014. Therefore, a reduction of the DD outbreaks in the northcentral area of the Iberian Peninsula is found during the period studied.
3. Overall, the mean DD contribution to AOD 440nm is 0.015 or 11.5%, while for the surface concentration PM 10 ,[5][6][7][8][9][10] this is 1.3 µg m -3 (11.8%), 0.55 µg m -3 (8.5) and 0.79 µg m -3 (16.1%), respectively. 4. The annual cycle of the DD contribution to aerosol load peaks in March, decreases in April-May, notably increases during summer months (the AOD curve has a local minimum in July), and experiences a progressive decline after summer (with a significant fall in September for the PM 10 curve) towards minimum values in winter.
The maximum DD contribution to AOD occurs in June and August close to 0.03, while the PM 10 maximum DD 10 contribution reaches ~2.4 µg m -3 in August. With careful inspection of all the information, the inventory can be a useful tool to develop and validate automated methodologies which use other instruments such as Raman lidars and ceilometers or use model forecasts. The comparison between different methodologies allows a more reliable estimation of uncertainties in DD detection and its contribution to 30 total aerosol load. Future studies based on this inventory will be focused on a global characterization of microphysical and radiative properties of desert dust including the evaluation of its radiative forcing over the study region. Therefore, these results are useful for assessing regional climate change studies linked to atmospheric aerosols because of the excellent clean background conditions of the area, which may be considered as one of the few sites/areas in south-western Europe with these conditions.

Acknowledgements
The authors are grateful to Spanish MINECO for the financial support of the FPI grant BES-2012-051868 and project 5 CGL2012-33576. Thanks are due to EMEP (especially to MAGRAMA and AEMET) and AERONET-PHOTONS-RIMA staff for providing observations and for the maintenance of the networks. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007(FP7/ -2013 under grant agreement Nr. 262254 [ACTRIS 2]. We also thank "Consejería de Fomento y Medio Ambiente" for their support to desert dust studies in "Castilla y León" region, as well as "Consejería de Educación of Junta de Castilla y León" for financing the project (VA100U14).