Trends in MODIS and AERONET derived aerosol optical thickness over Northern Europe

Abstract Long-term Aqua and Terra MODIS (MODerate resolution Imaging Spectroradiometer) Collections 5.1 and 6.1 (c051 and c061, respectively) aerosol data have been combined with AERONET (AERosol RObotic NETwork) ground-based sun photometer observations to examine trends in aerosol optical thickness (AOT, at 550 nm) over Northern Europe for the months April to September. For the 1927 and 1559 daily coincident measurements that were obtained for c051 and c061, respectively, MODIS AOT varied by 86 and 90%, respectively, within the predicted uncertainty of one standard deviation of the retrieval over land (ΔAOT = ±0.05 ± 0.15·AOT). For the coastal AERONET site Gustav Dalen Tower (GDT), Sweden, larger deviations were found for MODIS c051 and c061 (79% and 75%, respectively, within predicted uncertainty). The Baltic Sea provides substantially better statistical representation of AOT than the surrounding land areas and therefore favours the investigations of trends in AOT over the region. Negative trends of 1.5% and 1.2% per year in AOT, based on daily averaging, were found for the southwestern Baltic Sea from MODIS c051 and c061, respectively. This is in line with a decrease of 1.2% per year in AOT at the AERONET station Hamburg. For the western Gotland Basin area, Sweden, negative trends of 1.5%, 1.1% and 1.6% per year in AOT have been found for MODIS c051, MODIS c061 and AERONET GDT, respectively. The strongest trend of –1.8% per year in AOT was found for AERONET Belsk, Poland, which can be compared to –1.5% per day obtained from MODIS c051 over central Poland. The trends in MODIS and AERONET AOT are nearly all statistically significant at the 95% confidence level. The strongest aerosol sources are suggested to be located southwest, south and southeast of the investigation area, although the highest prevalence of pollution events is associated with air mass transport from southwest.


Introduction
Global warming is primarily a problem of too much carbon dioxide (CO 2 ) in the atmosphere. However, anthropogenic emissions of aerosols impact the Earth's radiation balance and climate as well (IPCC, 2013). The overall impacts from aerosols is a cooling of the Earth system. CO 2 is relatively stable in the atmosphere (with a residence time of about 100 years), and therefore evenly distributed over remote areas of the Earth.
Anthropogenic aerosols feature substantially shorter residence times in the lower atmosphere ($1 week) and cause mainly local or regional effects. At present day, the full impact of increased greenhouse gas concentrations on global temperature is not known as man-made aerosols mask the heating of the planet to an uncertain degree. Although the direct and indirect aerosol effects (in the cloud-free and cloudy atmosphere, respectively) have been known for decades, the estimates of the impacts of aerosols on the radiation balance are still associated with large uncertainties (IPCC, 2013). Even so, second to greenhouse gases, man-made aerosols have caused the largest anthropogenic forcing during the industrial era (Boucher et al., 2013). Characterizing global aerosol distributions and their changes over time is necessary for understanding present climate conditions and future climate change (e.g. IPCC, 2013), particularly as even natural aerosols have a large contribution to uncertainty in indirect forcing (Carslaw et al., 2013).
Anthropogenic emissions of aerosols and precursors have increased significantly since pre-industrial times (Lamarque et al., 2010;Granier et al., 2011;Smith et al., 2011). However, air quality mitigation measures in Europe and North America from the 1980s onwards have led to regional reduction in anthropogenic emissions; e.g. SO 2 emissions in Europe have declined with 73% between 198073% between and 200473% between (Vestreng et al., 2007Hand et al., 2012). Granier et al. (2011) report that anthropogenic emissions of NO x , CO and black carbon (BC) have decreased in Europe between 1980 and 2010 by 30%, 58% and 55%, respectively. From the Long-range Transmission of Air Pollutants inventory a decrease in SO 2 emissions for Europe from 21.65 Tg in 2000 to 16.68 Tg in 2010 has been estimated, which means a reduction of 23% (de Meij et al., 2012). The reduction in SO 2 for Europe given by EMEP for the period 1990-2010 is as much as 64% . For BC emissions the IPCC inventory shows a reduction of 11%  in Europe between 1900. However, between 2000 and 2010 a small increase of 3% (0.028 Tg) is found for this species. In the same report, the NOx emissions are about 10% lower  in 2010 compared to 2000. This reduction is suggested to be due to improved technology applied in road transport (e.g. catalytic converters and improved energy efficiency). For NH 3 emissions, the EMEP inventory shows a decrease with 2.5% between 2000 and 2010. This reduction is probably the outcome of the actions to combat eutrophication problems in northwestern Europe (de Meij et al., 2012). Furthermore, general there is a tendency that PM 2.5 levels measured in Europe have decreased over the period 2000-2014 (Guerreiro et al., 2016).
NASA has developed a suite of satellites known as the Earth Observation System to monitor important climate drivers, including the characterization of global aerosol properties and their changes over time. Two of these satellites are the twin MODerate resolution Imaging Spectroradiometer (MODIS, Salomonson et al., 1989) that provide measurements of the global distribution of aerosols from two polar orbiting satellites: Terra since 2000 and Aqua since 2002 (Remer et al., 2008). The swath width of 2330 km enables daily near-global coverage with MODIS, allowing for near-real-time monitoring of aerosol loading in the form of aerosol optical thickness (AOT) (Koren and Kaufman, 2004;Al-Saadi et al., 2005). The length of the MODIS AOT time series enables the investigation of temporal trends in aerosol burden. For example, Koukouli et al. (2010) found a negative trend in MODIS Aqua and Terra (MODIS A and MODIS T , respectively) Collection 5 AOT, with respect to the period 2000-2007, for the Southern Balkan/Eastern Mediterranean region. The latitudinal distribution of relative AOT tendencies derived from Level 3 MODIS global monthly AOT for the period 2000-2006 also shows negative trends within the range from -5.7% to about -3% per year at mid and high latitudes in the northern hemisphere (Kishcha et al., 2007). Based on data from both monthly MODIS Level 2 and 3 products (Collection 5), and from Aerosol Robotic Network (AERONET, Holben et al., 1998), de Meij et al. (2012 found significant negative trends in AOT for Europe for the period 2001-2009. Finally, sun-photometer data from 20 AERONET sites in Europe also show a decline in AOT: 11% over the period (Turnock et al., 2015. At the same time, an artificial tendency in the MODIS T Collection 5 AOT Level 2 product has been found over land, caused by a drift in the radiometric calibration (Levy et al., 2010;Zhang and Reid, 2010;Levy et al., 2013). Additionally, for the ocean Collection 5 retrievals a significant offset between MODIS A and MODIS T AOT has been found by Remer et al. (2008). Levy et al. (2013) therefore introduced calibration changes in Collection 6 due to the biases found with Collection 5. Even so, Levy et al. (2013) were not able to assess how the revised calibration will impact global trends and divergence of AOT derived from MODIS Terra and Aqua. Levy et al. (2010) state that global trends are difficult to quantify, due to complicated sampling patterns for different aerosol types and surface conditions. On the other hand, it is plausible that regional trends can be compared with ground-based or other satellites measurements (e.g. Papadimas et al., 2008;Karnieli et al., 2009).
In the present study, 13 and 15 years (2003-2015 and 2003-2017, respectively) of MODIS level 2 AOT data, derived with Collection 5.1 (c051) and Collection 6.1 (c061), respectively, have been used to estimate trends from daily and monthly (April-September) AOT for Northern Europe. For MODIS T also the years 2000, 2001 and 2002 are included in the investigation. The satellite derived time series of AOT have been compared to trends in AOT derived from AERONET measurements in the same region. Tesche et al. (2016) show that cloudfree conditions, required for passive remote sensing, only make up 20-40% of the daytime observations for the months from April to September in eastern middle Sweden. Therefore, we have investigated if the AOT availability is higher over the Baltic Sea compared to the surrounded land areas and what this means for the representativity of the data. Additionally, days with available MODIS AOTs obtained with c051 for Northern Europe have been compared to results obtained with c061. The availability in AOTs derived from one of the satellite platforms Aqua and Terra has as well been compared to AOTs derived from MODIS A and MODIS T observation combined over the investigation area. Collocated comparisons with AERONET have been performed to validate MODIS AOT. Finally, MODIS AOT has been combined with back trajectories calculations to investigate transport of polluted and clean air masses. The aim with the latter is to indirectly identify locations of aerosol sources that may have caused the changes found in AOT over Northern Europe.

MODIS aerosol retrievals
The MODIS Aqua and Terra Collection c051/c061 Level 2 standard 10-km products have been used here for best quality retrievals (quality flag ¼ 3) of AOT over land and ocean. Data were obtained through NASA Goddard Space Flight Center's Atmosphere Archive and Distribution System (http://ladsweb.nascom.nasa.gov). Detailed descriptions of the MODIS dark target algorithms for retrievals of AOT over land and ocean can be found in Remer et al. (2005), Levy et al. (2007), Levy et al. (2010) and Levy et al. (2013). After corrections for water vapor, ozone and CO 2 have been applied, the next step in the algorithm is to organize the reflectance at six wavelengths into 100 km 2 boxes of 20 Â 20 pixels at 500 m horizontal resolution. The ocean algorithm requires that all 400 pixels in the box are identified as ocean pixels. If any land is encountered the entire box is handled by the land algorithm. A test is used to filter out shallow-water pixels contaminated by underwater sediment (Li et al., 2003). Other tests filter out ice/snow pixels (Li et al., 2005), bright land scenes, and sun glint over water. The cloud screening is performed individually for the 400 pixels (Gao et al., 2002;Martins et al., 2002). Pixels that remain in a 10 km 2 box after the screening procedure are sorted according to their 858-nm brightness. The brightest and darkest 25% of the pixels over ocean are discarded, thereby leaving 50% of the cloud-free data for the retrieval. Over land, the brightest 50% and darkest 20% of the pixels are discarded, which means that AOT is retrieved if at least 10 of the 400 pixels in the original box remain after masking and filtering. The remaining pixels are assumed to represent conditions suitable for aerosol retrievals (Remer et al., 2012). The 1354 pixel wide MODIS granule, in fact, represents a 2330 km-wide swath at the surface. The tenth granule has an output size of 135 by 204 pixels, and the pixel size increases towards the swath edges. The MODIS land and ocean retrievals give AOT at 550 nm with expected error envelopes of DAOT ¼ ±0.05 ± 0.15 Ã AOT (Levy et al., 2010) and DAOT ¼ þ0.04 þ 0.1 Ã AOT, DAOT ¼ À0.02 À 0.1 Ã AOT (Levy et al., 2013), respectively, which arise from combined errors in assumed boundary conditions (e.g. surface reflectance, instrument calibration) and type of aerosol model (such as in single scattering albedo). In validation of MODIS against AERONET the expected error envelopes containing at least 67% (approximately one standard deviation) of the collocated AOT matches (Levy et al., 2013). The local equatorial crossing times are approximately 10:30 a.m. (descending node) and 01:30 p. m. (ascending node) local time for Terra and Aqua, respectively.
The aim in developing MODIS c061 was to improve the MODIS dark target aerosol product without a complete overhaul, thus, the basic theory, science and logic of the algorithm remain similar to c051 (Levy et al., 2013). However, coding bugs have been fixed, assumptions in the c051 algorithm have been reconsidered and diagnostic information has been added. As for the previous transition of going from Collection 4 to c051, the use of MODIS AOT data for determining aerosol climatology and trends impacts due to upstream calibration was also needed to be quantified (Levy et al., 2013). Although the aerosol algorithm was to remain unchanged from c051 to c061, the global aerosol product would be different, since the inputs are different. The new c061 product not only represents an update to the aerosol algorithm but also an update to all MODIS algorithms, including the calibration and cloud masking algorithms that produce the inputs to the aerosol algorithm (Levy et al., 2013). Note that MODIS AOT Level 2 standard 10-km products data derived with the c051 algorithm are only available up to 2015. Figure 1 shows the investigation areas that have been selected here to investigate AOT over Northern Europe. All MODIS Aqua and Terra AOTs within the large white dashed box in the figure, and that correspond to the months April to September of the periods listed in Table 2 are included in the investigation. To combine AOTs from MODIS Aqua and Terra we have introduced a common grid that consists of 0.25 latitude * 0.25 longitude pixels within the investigation areas. The three smaller areas denoted with white solid boxes in Fig. 1 have been selected to investigate trends in MODIS AOT. The marine areas of these boxes have been selected as there are more AOTs available over the Baltic Sea than over land (Section 4.3) and the positions are directly or nearly connected to the current ground-based AERONET stations. The European Centre for Medium Range Weather Forecasts (ECMWF) land-sea mask was used to exclude land pixels in the calculation of MODIS median AOT.

AOT from sun photometer measurements
MODIS AOT retrieved over land has been compared to AERONET (version 3) level 2.0 data (quality assured) from Hamburg (53.6 N,10.0 E),Germany,and Belsk (51.8 N,20.8 E), Poland, for the periods 2003-2016 and 2002-2016 (not 2003), respectively. In addition, over the ocean MODIS AOT has been compared to AERONET (version 3) level 2.0 data from Gustav Dalen Tower (GDT; 58.6 N, 17.5 E), Sweden, for the period 2005-2017. For GDT data for April, May, August and September are not available for all years investigated, while for Belsk the data set covers only the period from April to July in 2012.
Information about the CIMEL sun photometers operated at these sites can be found at http://aeronet.gsfc. nasa.gov. AERONET data used for this study include AOT at 500 nm as well as the Ångstr€ om exponent a (440/ 675 nm). These values were recorded every 15 minutes and automatically cloud screened (Smirnov et al., 2000). AERONET-derived estimates of spectral AOT are expected to be accurate within < ±0.01 for wavelengths larger than 440 nm (e.g. Holben et al., 1998).
Sun photometer data from the considered sites have been used as well to investigate trends in AOT, but also to validate MODIS AOT. For Belsk and Hamburg the Ångstr€ om power law was used to convert AERONET   Table 1 for additional information.
AOT at 500 nm to the wavelength of 550 nm, for which AOT is retrieved from MODIS observations. The AOT at GDT is retrieved directly at 550 nm. Nine satellite pixels around each site were used in the validation of MODIS AOT, although only cases with an availability of at least five AOT values were included in the comparison. The ground-based sun photometer measurements have been averaged over 2.5 hours around the time of each satellite overpass. The relative (normalized) root mean square deviation (RMSD rel ) was determined for collocated satellite and ground-based AOT values, averaged according to the procedure described in (Mishchenko et al., 2010).

HYSPLIT trajectory model
Hourly 5-days back trajectories were calculated using the HYSPLIT model (Draxler and Hess, 2003) for April to September of the period 2003-2017. The back trajectories start at 500 m asl from western Gotland Basin (wGB,58.38 N,17.88 E) and southwestern Baltic Sea (swBS,54.75 N,13.25 E). The trajectory calculations are based on meteorological data from the NCEP/NCAR reanalysis at 2.5 resolution (http://ready.arl.noaa.gov/archives.php).
Complementary information about the current coordinates of the investigations areas, as well as of the AERONET and HYSPLIT sites, investigation periods and a summary of the present objectives are found in Table 1.

Method to investigate trends in AOT
Trends in AOT have been estimated with respect to MODIS A and MODIS T observations separated and combined over the investigation area. A screening procedure has been introduced here for the MODIS observations to assure a sufficient number of AOT values are available for representing a full day in the estimations of trends. Days with the number of AOT values equal and lower than 5% (5th percentile) of all available AOTs within the investigation area and corresponding to the investigation Temporal autocorrelation can be a confounding factor when analysing trends in time series data (Weatherhead et al., 1998;von Storch and Zwiers, 1999). Temporal autocorrelation, if present, would effectively reduce the number of independent observations in a given time series. Autocorrelation has been taken in consideration in the trend analyses by using a stable seasonal adjustment on the current time series of AOT (Brockwell and Davis, 2002). In addition, temporal autocorrelation has been taken in consideration as well by calculating the autocorrelation function (Box et al., 1994) for all cases investigated in the present study. The method described in von Storch and Zwiers (1999) was used to calculate adjusted values for the degrees of freedom in our samples. Using this method, the adjusted degrees of freedom were substantially lower than the number of samples, particularly for the trends estimated based on daily AOTs: more than one order in magnitude. The significance of the trends were tested based on a z-test and the adjusted degrees of freedom values.
Some degree of correlation, R 2 ¼ 0.31 and R 2 ¼ 0.16, is found in the time series of MODIS c061 daily AOT between swBS and wGB as well as between swBS and cP, respectively. Simultaneously, significant differences in the availability in AOT between the investigation areas occur: 79% and 52% daily coverage for wGB and cP, respectively, with respect to swBS, This means that dealing with spatial autocorrelation is not fully relevant here due to widely differing data coverage and has therefore not been taken in consideration in the present investigation.

Method to investigate spatial transport pattern of polluted and clean air masses
To investigate spatial patterns of air mass transport we used 5-day back-trajectories with daily starts at 09:00, 10:00, 11:00, 12:00 and 13:00 UTC, for April to September of the period 2003-2017. The coordinates of the back-trajectory hourly positions, were projected to a polar grid with the origin located at the trajectory receptor location within swBS box (Fig. 1). The probability of occurrence is defined as the number of trajectories that have passed a cell divided by the total number of trajectories used. In the case, that a trajectory crossed a grid cell more than once only one hit was registered. The fraction of hits of the days corresponding the lowest percentile (25th) and highest percentile (75th) of the MODIS c061 AOTs (AOT 25 ¼ 0.12 and AOT 75 ¼ 0.25, respectively) were compared in order to investigate the spatial transport patterns of cleaner and more polluted air masses. The corresponding areas for MODIS daily AOT estimates are swBS (Fig. 1).
Furthermore, the contribution of each of 16 azimuthal sectors to the total long-term AOT at wGB and swBS was calculated. This was done by linking MODIS c061 daily AOT with the fraction of time that each of 5-day back-trajectories of that day was present in each sector. Only trajectory heights that were within the mixed layer were accounted for, since the source of the particulate matter that affects AOT is assumed to be at the surface. The mean AOT for each sector, weighted by the residence time of the trajectory points in each sector, was also calculated. This was done in order to reveal the directional dependence of AOT in each of the two investigation sites.  Figure 2 shows that MODIS AOT, for both c051 and c061, are well within the expected error envelopes when compared to AOT derived from sun-photometer measurements at Hamburg and Belsk. The figure also shows that as many as 24% more matches are available for c051 compared to c061 even though the latter collection also includes data from 2016. When comparing c051 AOT against AERONET, 89% of 901 collocated matches are found to be within the predicted uncertainties for MODIS A , while there are only 83% for MODIS T (1026 collocated matches). For c061 applied on MODIS A and MODIS T observations, 92% (687 collocations) and 88% (872 collocations) of matches, respectively, are found to be within the predicted uncertainties. Figure 3 shows a larger deviation in AOT between MODIS and the sun photometer at GDT compared to the validation over land. When comparing c051 AOT against AERONET for MODIS A and MODIS T separated 80% and 78% of matches (382 and 505 collocations, respectively), respectively, are found to be within the predicted uncertainties. For c061, 80% and 70% of matches are found to be within the predicted uncertainties for MODIS A and MODIS T , respectively (728 and 681 collocations, respectively). Furthermore, Fig. 3 shows a factor of 1.6 higher 6 P. GLANTZ ET AL.

Comparison of MODIS AOT against AERONET measurements
number of collected points for c061 than c051 for this coastal site. In Section 3.2 we will show the advantage of using measurements from both MODIS A and MODIS T to estimate AOT over the investigation area. Figure 4 shows relative frequency histograms of AERONET AOT and cumulative distribution curves of MODIS and AERONET AOT, subdivided according to the ground-based stations Belsk/Hamburg and GDT as well as for c051 and c061. The cumulative distribution curves shown in the figures have been subdivided according to MODIS A and MODIS T . For Belsk/Hamburg the difference in median AOT of the cumulative distribution between satellite and ground-based remote sensing is relatively small, although larger deviation is found for c061 applied on MODIS T . For GDT larger deviation between MODIS and AERONET is found.  3.2. Influences of cloud formation and sun glint effects on the availability in AOTs Figure 5(a, c) shows MODIS Terra and Aqua AOT c061 scenes, respectively, over Northern Europe for 2 July 2014. Figure 5(b, d) shows the corresponding visible composite images, in which cloudy and cloud-free areas can be identified. AOTs derived from MODIS Terra observations over the North Sea around 10:30 UTC are denoted with the red dotted oval in Fig. 5(a). For the same area Fig. 5(c) shows that MODIS A that passed this area around 12:20 UTC in the same day has not produced any AOTs. This is likely due to sun glint effects, specular reflection of light from water surfaces, which means that AOTs from MODIS A were screened out. This problem arises when sunlight is reflected off the water surface at the same angle that the satellite sensor views. For MODIS A this appears along and near the nadir view (left hand side of the satellite track). The orange solid oval suggests that the opposite occurs over the Baltic Sea: sun glint in MODIS Terra observations occurs on the right hand side of the nadir view. Thus, using combined observations from both MODIS A and MODIS T over the current marine area means that the availability in daily AOTs increases by as much as40% compared to using data from just one of the two platforms (Table 3). Figure 6 visualizes the effect of formation of convective clouds over land on the AOT retrieval. The MODIS A overpass at 1235 UTC shows increased cloud coverage in the marked area (orange solid oval) compared with the nearly cloud-free MODIS T scene from 09:10 UTC earlier that day. For such condition, it can thus be expected that MODIS T associated with earlier overpasses over the current investigation area lead to a better AOT coverage than for MODIS A . show that the major part of the investigation area is associated with higher availability in AOTs for the combined data set. Even so, the figure shows that very low number of days with available AOTs appears over Sweden, particularly for MODIS A observations in April. The latter may be due to the development of shallow convective clouds between the overpasses by Terra and Aqua. Nevertheless, along parts of the coastal land higher number of days with available AOTs is found. Better statistics are suggested also for two belts in the east-west direction located in the middle part of eastern Sweden, which is most clearly shown in Fig. 7(f). Furthermore, Fig. 7 shows that substantially higher numbers of days with available AOTs are found over water rather than over land. This is not valid, however, along the coasts and for the Danish islands. When comparing the findings in Fig. 7 with results obtained with MODIS c061 large differences in AOT are found between the two collections for the retrievals over land. Figure 8(a , b) shows that the majority of the land areas are associated with significantly higher number of daily AOT for c051 compared to c061. However, along a large part of the coastline (including over GDT), over the large Swedish inland lakes and between the Danish islands the opposite appears. In addition, large difference in available AOTs appears when comparing MODIS c051 and c061 scenes corresponding to 22 September 2006 (Figs. 1 and 9(a), respectively). By comparing these results with MOD03 c061 land sea mask ( Fig. 9(c)), e.g. for Southern Sweden and Northern Germany, it is obvious that this updated screening plays a role for the difference found in the statistical representation of AOT. The MODIS visible composite pictures in Fig. 9(c, d) suggests cloud free conditions over these areas.

Availability in AOTs derived from MODIS aqua and terra observations
Since the majority of the aerosols that contribute to AOT in the investigation area are probably accumulation-mode particles (Levy et al., 2010), relatively small spatial variations are expected for Northern Europe (e.g. Tesche et al., 2016). This, together with the finding that Baltic Sea provides better statistical representation of AOT than the surrounding land, support a focus on specific regions over the Baltic Sea in order to estimate trends in AOT (Section 3.5) as well as to investigate transport patterns of clean and polluted air masses (Section 3.6).  Figure 10 shows MODIS c051 Aqua/Terra monthly averaged AOT fields over the investigation area, for the period 2003-2015. Note that the number of days used for this averaging varies considerably between the different regions, as indicated in Fig. 7. The representativeness of the mean AOD over land in the northern part of the area of interest may therefore be limited. A better balance of the number of days with available AOT over land and water for July compared to the other months (not shown) seems to produce more coherent AOT fields at lower and higher latitudes ( Fig.  10(d)). Despite differences in AOT availability, Fig. 10 suggests a north-south gradient in AOT. The latter is presumably because the abundance of aerosol sources are located south of the present investigation area.

MODIS c051 monthly AOT fields
Mean AOT values that are substantially higher than a typical value corresponding to a background aerosol ($0.1) are also found, particularly in the southern part of the investigation area. However, the figures indicate that anthropogenic aerosols seem to influence the investigation area to a lesser degree in September compared to the other months investigated. Lower humidity and consequently lower hygroscopic growth of the aerosols in fall may also play a role.

Trends in satellite and ground-based derived AOT
Time series of monthly median AOT over the periods 2003-2015 and 2003-2017, obtained with c051 and c061, respectively, applied on MODIS A/T observations over the swBS as well as AERONET measurements at Hamburg, are shown in Fig. 11. For MODIS T trends in AOT have been estimated also based on the years 2000, 2001 and 2002. The solid and dashed lines denote linear fits and corresponding 95% confidence intervals for the slope, respectively. Trend values (in % per year) estimated for monthly AOT with respect to the present investigation period are also shown. Similarly, Figs. 12 and 13 show the AOT trends for central Poland (cP) and wGB with the sun-photometer measurements at Belsk and GDT, respectively.
In addition, Table 2 presents a summary of statistics obtained from the results shown in Figs. 11-13, with additional information on the separate analysis of MODIS A and MODIS T as well as on trends derived from daily AOT values. Regression statistics based on daily medians instead of monthly median values means higher degrees of freedom. Since bias caused by calibration degradation has been found for MODIS Terra c051 (Section 4.2) no values of the slope and y-intercept are presented in Table 2 for this platform. Table 2 shows that negative AOT trends are found for all cases and the majority of these are statistically significant at the 95% confidence level. The trends estimated over the three investigation areas are in line with negative trends found at the AERONET stations, although the differences between the MODIS trends estimated for cP area are very large. Note that GDT is associated with substantially lower availability in daily AOT in April, May, August and September due to missing data in some years. The data-coverage threshold included in this study was selected to avoid too low numbers of AOT values for representing a day (Section 2.4). This somewhat reduced the overall availability in daily AOT values. By including and excluding the threshold for all MODIS cases (Table 2) gives RMSD rel and RMSD of 13% and 0.16% per year, respectively. Thus, the difference is relatively small, but not negligible.
A summary of the data availability for each month and dataset is given in Table 3. Tables 2 and 3 confirm that the number of days with available AOTs are on the whole higher for c051 than c061, and it is particularly large for retrievals of AOT over land.

Air mass transport and aerosol sources
A small ensemble of 5-day back-trajectory was calculated for each day with available MODIS A/T c061 AOTs for the swBS region. The receptor point is indicated in Fig. 1. Figure 14 displays a probability map of the trajectories crossing the grid cells during the period of interest. Panel a) comprises only trajectories from the days of the lower quartile of AOT (<0.11), while panel b) is based on the highest quartile (AOT >0.25). This provides an insight into the difference in the general transport patterns between the days with the lowest and highest AOT. Both panels imply that westerly and north westerly winds are dominant over the investigation area, but for cleaner days, the air masses were mainly transported over the North Atlantic and Scandinavia before reaching swBS. The trajectories associated with higher AOT have a significant southwestern component. This is an indication that anthropogenic aerosol sources located in Western Europe, including the English Channel, likely influence the swBS. The higher probability of occurrence west of Great Britain does not imply that the Atlantic is an important contributor to the swBS AOT. It is more likely an indication of the position of the trajectory points before crossing the source regions. Figure 15 shows the directional dependence of the mean MODIS AOT for swBS and wGB. It is based on the residence time of the trajectories in each sector and the MODIS c061 AOT for both Aqua and Terra. The colour-shading indicates that the southerly sectors for both locations are associated with the highest mean AOT, which is around double of the mean AOT associated with the northern sectors. It is also clear that the mean AOT is higher in all corresponding swBS versus wGB sectors. This is likely due to swBS being closer the continental Europe and busy shipping lanes (Aulinger et al., 2016), where air masses with higher AOT arrive from. Furthermore, the length of the sectors in Fig. 15  indicates their contributions to the total long-term AOT at both locations. Despite not being associated with the highest mean AOT, the high frequency of air-masses arriving from the two sectors around the western line dominates, and yields a contribution of 25%-30% to the total AOT in a long-term perspective. The north-eastern sectors contribute the least to the total AOT because trajectories arriving from these sectors are quite infrequent, and they are also associated with relatively low AOT loading.

Validation and availability of MODIS AOT
AOT derived with c051 and c061 from MODIS observations over Northern Europe, for April-September, has been validated against AOT obtained from sun photometer measurements at two land sites and one coastal station. MODIS A/T 550-nm AOTs derived with c051 vary within the expected uncertainties of the MODIS retrievals over land (DAOT ¼ ±0.05 ± 0.15 Ã AOT) for 86% of the  1927 compared cases, although the best agreement occurs for MODIS Aqua (89%). Values of R 2 ¼ 0.86 and RMSD rel ¼ 31% were found for the MODIS c051 and ground-based collocated cases. The present results are in agreement with Glantz and Tesche (2012) who included 20 AERONET stations in Europe in the validation of MODIS c051 and Tesche et al. (2016) who focussed on long-term MODIS observations over the Stockholm region. In a global validation of the MODIS c051 darktarget aerosol products against AERONET, Levy et al. (2010) conclude that MODIS compares best (negligible bias and good correlation), for light aerosol loading (AOT <0.15), over sites that are both moderately 'dark' and moderately 'green'. Such generally vegetated sites occur over the Eastern United States and Southern Africa, but also over Western Europe. For heavy aerosol loading (AOT >0.4) the best agreement occurs among others for Western Europe, where the fine-mode aerosol dominate (Levy et al., 2010). MODIS A/T 550-nm AOTs derived with c061 was found to vary within the expected uncertainties of the MODIS retrievals over land for 90% of the compared cases. Values of R 2 ¼ 0.89 and RMSD rel ¼ 27% were found for this collection. However, MODIS T 550-nm median AOT derived with c061 over land was found to be 15% higher than AERONET median AOT. In addition, it was found here that collocated matches over land between MODIS and AERONET are substantially higher for c051 than c061. Using c061 and MODIS 3-km in the study by Tesche et al. (2016) also reduces the available AOTs over Sweden compared to MODIS c051. The MODIS c061 retrieval uses more aggressive tests for spatial variability of reflectance and normalized difference vegetation index (NDVI) to mask clouds and inland water bodies (Levy et al., 2013). Livingston et al. (2014) recommend exploring whether different thresholds on cloud and NDVI masks can expand the MODIS aerosol data set in condition of smoke plumes, without a risk for contamination by clouds and highly reflecting surfaces. The present study suggest that the consequences of this new NDVI mask introduced should be further investigated also for situations with lower aerosol loadings. MODIS A/T 550-nm AOTs derived with c051 and c061 were found to vary within the expected uncertainties of the MODIS retrievals (DAOT ¼ þ0.04 þ 0.1 AOT,   14 P. GLANTZ ET AL. DAOT ¼ -0.02 -0.1 AOT) over the coastal site GDT in the wGB, Sweden, for 79% and 75%, respectively, of the compared cases. A reason for the lower agreement in AOT between MODIS and this AERONET station is probably due to the location of the GDT (M elin, 2011; Levy et al., 2013;M elin et al., 2013): $30 km from the east coast of Sweden, with a minor archipelago in between. By combining AOT derived from both MODIS A and MODIS T the present results show that the availability in daily AOT increases substantially compared to using data from only one of the platforms. The present results also show that the number of days with available AOT is substantially higher over the Baltic Sea (away from the coastlines) compared to the nearby land areas. The latter is confirmed in Table 3 by comparing mean number of days with available AOTs between the AERONET stations GDT (ocean) and Hamburg/Belsk (land) for June and July. Higher availability in AOTs over the Baltic Sea is likely due to fewer clouds formed over water compared to the surrounding land.
In the present study it was found however that lower availability of MODIS AOTs occurs along many parts of the coastlines and between islands (e.g. in Denmark). An improper representation of the coastline by MODIS c051 has also been found in the study by Tesche et al. (2016). Higher availability of MODIS AOT over the present areas was in any case found for c061 than c051. Even so, M elin et al. (2013) suggest that further improvements of the atmospheric correction in coastal waters may require an additional level of detail to describe the complex aerosol mixtures more accurately, and specific water biooptical models for an appropriate representation of the boundary condition. Furthermore, very low numbers of days with available AOTs were found over Sweden, particularly in April. Even so, along some of the coastlines higher numbers of days with available AOTs were found. Higher availability in AOT was also found for two belts in the east-west direction located in the middle part of eastern Sweden. Better AOT statistics for the northern belt is likely due to the location of the third largest lake in Sweden, called M€ alaren. Better AOT statistics corresponding to these two belts and along the land coast are confirmed by available sun hours over land (not shown). The latter is obtained based on pyrheliometer measurements combined with daily model calculation carried out by the Swedish Meteorological and Hydrological Institute.
When comparing number of days with available AOTs between MODIS c051 and c061, a substantially higher availability was found for retrievals over land with the former collection (Fig. 8). The difference exceed 15% in April and July for a large part of the land area south of the Baltic Sea, while it is even larger over Sweden. Caroll et al. (2017) shows that the updated land mask in c061 identifies a lot more smaller inland water bodies. This finding in combination with the more aggressive tests for spatial variability of reflectance and normalized difference vegetation index (NDVI) to mask clouds and inland water bodies (Levy et al., 2013), introduced in c061, likely explain this difference in AOT statistics. Simultaneously as this large difference in data coverage occurs between the two collections, the validation against AERONET shows only a slightly better agreement for MODIS c061. These findings suggest that the c061 algorithm needs to be further developed so that more of inland pixels adjacent to water bodies can be included in the retrievals of AOT. If accuracy in the retrievals of AOT for mixed pixels can be obtained we believe that passive remote sensing with MODIS can be even more useful to estimate longterm trends in AOT.

Trends in MODIS c051/c061 and AERONET AOT
Changes have been made with MODIS c061 considering the algorithm used in the retrievals of AOT, but also to address problems with calibration degradation in previous collections. A drift in c051 AOT was found for MODIS T land retrievals (Remer et al., 2008;Levy et al., 2010). Levy et al. (2010) validated this collection against AERONET AOT and found an overestimation of about 0.005 prior to 2004 and an underestimation of the same magnitude afterwards for MODIS T . This is likely caused by a drift in the radiometric calibration (Zhang and Reid, 2010;Levy et al., 2013). For the ocean, c051 retrievals a minor offset has also been found between MODIS A and Fig. 13. Same as Fig. 11, but for MODIS derived AOTs over wGB and the AERONET GDT site.

TRENDS IN MODIS AND AERONET
MODIS T AOT (about 0.015 higher for the latter platform, Remer et al., 2008). Thus, these biases associated with c051 applied on MODIS T AOT observations may, at least partly, explain the differences found here in AOT between the two platforms and collections. This suggest that results obtained with MODS c061 should be given priority. However, since the present study also found large differences in the availability in AOT between the two collections trends in AOT obtained with MODIS c051 are presented as well.
Trends in MODIS and AERONET AOT have been investigated with respect to the months April-September.
Since better availability in daily AOTs occurs over the Baltic Sea compared to the surrounded land areas this finding has been used to estimate changes in MODIS AOT for the period 2000-2017. Negative trends within the ranges 0.84%-1.70% and 1.10%-1.40% per year in AOT based on daily values were found for the swBS with MODIS c051 and c061, respectively. This can be compared to a decrease of 1.50% per year in AOT at the AERONET Hamburg station. On the whole stronger negative trends in AOT are found over the cP area and wGB for MODIS c051 and AERONET. However, the absolute values of AOT are substantially lower for  Table 1. AERONET GDT compared to MODIS AOT, which may be due to missing data in some years at this station, particularly for April. Particularly over cP, weaker negative trends in AOT are obtained with MODIS Collection 6.1. The negative trends estimated with MODIS c051, c061 and AERONET are nearly all statistically significant at the 95% confidence level. The present decrease in AOT during the current investigation period results in a reduction of the direct aerosol cooling effect during the high-insolation period in Northern Europe. The representation of MODIS AOT data are poor over Northern Europe during the low-light months from October to March, which means that these month are not included in the present investigation. Thus, this induces uncertainties in the present trend estimates. Even so, in subsequent investigations we plan to estimate if the decrease in AOT here has contributed to the "extra" warming that is observed over northern Europe. For the swBS and wGB areas ECMWF mean surface temperature has increased by as much as 1.8 and 1.9 , respectively, over the period 1979-2017 (April-September). The negative trends estimated in the present study are in line with a decline in AOT with 1.1% per year that has been obtained from observations at 20 AERONET stations in Europe within the period 2000-2009 (Turnock et al., 2015). Substantially stronger trends in AOT (555 nm), -1.8%, -2.1% and -1.9% per year, are obtained from Multi-angle Imaging SpectroRadiometer (MISR) observations over the swBS, cP and wGB, respectively, with respect to the period 2000-2016. However, the daily data coverage is very low for MISR due to the narrow swath width. For example, for the swBS and period from April 2000 to May 2017 (April to September) this means an availability in MSR daily AOT as low as 12%. The MISR AOT can be accessed through Giovanni (https:// giovanni.gsfc.nasa.gov/giovanni/). Moreover, aerosols with a dry diameter <2.5 mm present in the boundary layer normally contribute significantly to the total attenuation. The total mass of these particles, PM2.5, is commonly measured for air-quality monitoring and can be accessed through EBAS (http:// ebas.nilu.no). At several Swedish sites (e.g. Vavihill, Aspvreten, Norr-Malma) we found that PM2.5 has been declining by 3%-5% per year over the last decade. The numbers slightly vary depending on the site and the recorded period. Similar declines are seen in rural sites in Northern Germany (Waldhof) and Poland (Diable Gora). This indicates that at least part of the present negative trend in AOT in Northern Europe is due to a cleaner boundary layer.
Elevated aerosol layers that may induce uncertainties in the present AOT trends have not been investigated in the present study. At least over middle Sweden, it has however been found from spaceborne CALIPSO lidar observations that elevated aerosol layers associated with long-range transport events are rarely present (Tesche et al., 2016). The same has been observed from aerosol lidar measurements in Stockholm, Sweden, which were available between 2011 and 2014 (Baars et al., 2016). Furthermore, the contribution from stratospheric AOT on the column AOT was low an increased during the period 2002-2011 (Glantz et al. 2014), which means that volcanic aerosols (Vernier et al., 2011) do not contribute to the negative trends found here.

Transport of polluted and clean air masses
Finally, by combining 5-days back-trajectories with MODIS A/T AOT derived with c061 transport of clean and polluted air masses have been investigated. The results show that swBS is associated with higher mean AOTs compared to wGB for all transport directions. For both investigation areas, the cleanest situations appear when air masses are transported over Scandinavia, including the surrounding seas. However, the higher AOTs found for swBS, compared to wGB, also for northerly winds, may be due to local ship emissions. The strongest aerosol sources are suggested to be located southwest, south and southeast of the investigation area, although the highest prevalence of these events is associated with the former direction.
To quantitatively identify aerosol sources is beyond the scope of this paper. Even so, beside Western Europe the present results suggest that aerosols from commercial ships cannot be excluded as an important source of aerosols for Northern Europe (Aulinger et al., 2016). The English Channel is one of the busiest shipping lanes in the world. Strong sources of aerosols, although not so frequently occurring, are also suggested to be located east of the investigation area. Moving the site for the backtrajectories to Minsk in Belarus, east of the investigation area, and combining it with the 75th percentile of AOT suggest a substantially higher frequency of transport of air masses with high aerosol content from east (not shown). This is likely explained by agricultural or forest fires in Ukraine or Russia (Stohl et al., 2007;Cook et al., 2008;Glantz et al., 2009). of ECMWF data sets. We wish to thank the anonymous reviewer for the helpful comments. Finally, the authors are very thankful for the use of in situ aerosol data that are measured by the Department of Analytical Chemistry and Environmental Science (ACES) at Stockholm University.

Disclosure statement
No potential conflict of interest was reported by the authors.