A novel method of identifying and analysing oil smoke plumes based on synergic satellite data

Black carbon aerosols are the second largest contributor to global warming while also being linked to respiratory and cardiovascular disease. These particles are generally found in smoke plumes originating from biomass burning and fossil fuel combustion. They are also heavily concentrated in smoke plumes originating from oil fires exhibiting the largest ratio of 15 black carbon to organic carbon. In this study, we identified and analyzed oil smoke plumes derived from 30 major industrial events within a 12-year timeframe. To our knowledge, this is the first study of its kind that utilized a synergetic approach based on satellite remote sensing techniques. One objective of this study is to highlight the importance of satellite remote sensing techniques in identifying these types of events. As opposed to ground stations, satellite data offers access to remote areas all over the globe which would otherwise be very difficult to reach. Satellite data offers access to these events which, 20 as seen in this study, are mainly located in war prone or hazardous areas. This study focuses on the use of MODIS (Moderate Resolution Imaging Spectroradiometer) and CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations) products regarding these types of aerosol while also highlighting their intrinsic limitations. By using data from both MODIS instruments onboard Terra and Aqua satellites we addressed the temporal evolution of the smoke plume while assessing Lidar specific properties and plume elevation using CALIPSO data. We present several aerosol properties in the form of 25 plume specific averaged values. The MODIS ocean algorithms were successful in retrieving aerosol properties which, on average, ranged from -0.06 to 0.16 for plume specific AOD, 0.18 to 1.25 for Ångström exponent and 0.29 to 1.73 μm for the effective radius. CALIPSO measurements showed values of plume AOD ranging from 0 to 0.14 (532 nm) and 0 to 0.13 (1064 nm) except for one event where AOD values showed 1.52 (532 nm) and 1.42 (1064 nm). AE values ranged from 0.11 to 0.33 which were in agreeance with MODIS values. A large discrepancy can be found in one event where CALIPSO 30 measured AOD values 5 times higher than MODIS. This event also produced the largest lidar ratio at 109 sr (532 nm) and 86 (1064 nm). Other lidar ratio values ranged from 37 to 55 sr however these unconstrained solutions were obtained for the entire layer of which the plumes were a part of and thus did not reflect specific plume conditions. Particulate backscatter https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c © Author(s) 2021. CC BY 4.0 License.

values ranged from 0.002 to 0.0017km -1 sr -1 while extinction coefficient values ranged from 0.10 to 1.65 km -1 . On average backscatter and extinction coefficient values were 2 to 9 times higher than local background values. Particulate 35 depolarization ratios ranged from 0.11 to 0.15 in 4 out of 6 cases while the remaining two ranged from 0.27 to 0.32 where dust was highly dominant. The values represented in this study are in good agreement with similar studies that used ground based and flight measurements. We believe that MODIS values are a conservative estimation of plume AOD since MODIS algorithms rely on general aerosol models and various atmospheric conditions within the look-up tables which do not reflect the highly absorbing nature of these smoke plumes. CALIPSO measurements are heavily dependent on lidar ratios which are 40 not directly measured if plumes within the planetary boundary layer. We also believe that AOD values based on CALIPSO measurements are conservative in nature since heavy absorbing smoke would yield larger lidar ratios and AOD values.
Based on this study we conclude that the MODIS land algorithms are not yet suited for retrieving aerosol properties for these types of smoke plumes due to the strong absorbing properties of these aerosols. We believe that these types of studies are a strong indicator for the need of improved aerosol models and retrieval algorithms. 45

Introduction
Atmospheric aerosols are chemically complex mixtures of solid and liquid particles dynamically suspended in air.
They originate from both natural and anthropogenic emissions. More common naturally occurring aerosols can be observed in the form of fog, dust, sea salt spray, biological exudates and grey smoke (biomass burning). Haze, smog and black smoke are typically a result of industrial and transportation activities (Stocker et al., 2014;Wei et al., 2020). Distinct species such as 50 black carbon (BC), organic carbon (OC), sulphates, nitrates, trace elements, sea salt, mineral dust and biological matter suffer atmospheric alteration resulting in different combinations of compounds. Defining aerosol types is a difficult task as they possess a large degree of variance in composition and concentration due to different atmospheric residence times, dry deposition and wet scavenging, emission rates and sources, transport trajectories and seasonal variability (Dutkiewicz et al., 2009;Li et al., 2015;Samset et al., 2018). 55 Health effects associated with both short and long-term exposure to aerosol have been widely documented in scientific literature (Brauer et al., 2015;Laumbach and Kipen, 2012;Zhang and Batterman, 2013;Pascal et al., 2013;Lee et al., 2014;Guarnieri and Balmes, 2014;Zhang et al., 2017). Aerosols have been linked to respiratory and cardiovascular diseases due to fine particulate matter (PM2.5; <2.5 μm in diameter) that can penetrate the lungs resulting in increased rates of morbidity and mortality (Brunekreef and Holgate, 2002;Pope III et al., 2002;Lim et al., 2012;Beelen et al., 2014;Hoek 60 et al., 2013).
Global circulation of aerosols is a known transport vector of minerals and nutrients to the biosphere (McTainsh and Strong, 2007;Maher et al., 2010;Lequy et al., 2012). Aerosols have a direct effect on radiation distribution by scattering, absorbing and emitting light thorough the atmosphere. In addition they can affect the climate system through indirect effects acting as cloud condensation nuclei, impacting cloud lifetime and properties, atmospheric stability and precipitation factors 65 https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License. (Draxler et al., 1994;Limaye et al., 1992). Several international research teams conducted extensive field campaigns 100 concentrating their efforts on atmospheric, environmental and health related issues focusing on the potential impacts on global climate (Sadiq and McCain, 1993;Husain, Tahir, 1995;World Meteorological Organization, 1993). Studies on health effects related to the Kuwait oil fires suggest that smoke exposure lead to acute respiratory illnesses with some suspecting long-term effects (Etzel and Ashley, 1994;Brain et al., 1998;Smith, 2002;Lange et al., 2002;Kelsall, 2004;Heller, 2011;Barth et al., 2016). Valuable atmospheric data was also collected from smaller events such as oil depot fires, most notably 105 the Bouncefield incident on 11 December 2005. A number of explosions led to a large fire engulfing 20 storage tanks until the 15 th of December. The fire burned 58000 tons of fuel while injecting a large smoke plume above the boundary layer at 3000 m (Vautard et al., 2007;Health and Safety Executive and Buncefield Major Incident Investigation Board (Great Britain), 2008). An initial report on air quality concluded that the smoke plume remained aloft over cold and stable atmospheric layers, thus reducing the potential impacts at ground level (Targa et al., 2006). Health studies related to the 110 event concluded no long-term impacts on people exposed to the smoke however acute respiratory symptoms were reported (Hoek et al., 2007;Morgan et al., 2008).

Objectives 125
One objective of this study is to highlight the importance of satellite remote sensing techniques in identifying these types of events. As opposed to ground-based data, satellite data offers access to remote areas all over the globe which would otherwise be very difficult to achieve. As seen in Table 1, oil installations may be situated in desert areas, at sea or in secluded locations far away from air quality monitoring stations or AERONET (AErosol RObotic NETwork) sites (Holben et al., 1998). In addition to this advantage, a synergistic approach using different types of satellite instruments can offer 130 three-dimensional space coverage. While in situ, ground stations and modelling tools are viable options for smoke plume research these methods have limitations, in areas prone to armed conflicts or posing high health risks. Out of the https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License. aforementioned events only the event in Kiev was analysed by different techniques, as seen in local AERONET data (section 3). It goes without saying that retrieving optical and microphysical properties of petrochemical burnings may be challenging in most cases even with this approach. This study will focus on the use of MODIS and CALIPSO aerosol products regarding 135 these types of aerosol while also highlighting their limitations. By using data from both MODIS instruments on-board Terra and Aqua satellites we addressed the temporal evolution of the smoke plume while assessing Lidar specific properties and plume elevation using CALIPSO data. The low number of studies on petrochemical smoke plumes, especially in the last decade, further encourages us to address these issues. While biomass burning and industrial haze are abundantly discussed in scientific literature, the same cannot be said for petrochemical smoke plumes resulting from major technological accidents. 140 To our knowledge, we have not identified any similar studies focused specifically on retrieving aerosol properties from major petrochemical accidents by using synergistic satellite techniques.

Event synopsis
This section summarizes a collection of events ranging from 2008 to 2019 that were successfully identified by satellite remote sensing techniques. Table 1 also provides the coordinates and the number of MODIS observation for each of 145 the events covered in this study.
Events similar in nature to the Kuwait oil fires took place in Northern Iraq as oil fields at Qayyara and Najma were intentionally set ablaze by Islamic State of Iraq and Syria (ISIS) militants in an attempt to deter coalition air strikes. The first fires were detected east of Baiji in early January 2016. Other oil wells burned intermittently form May to June close to Mosul and Kirkuk. The bulk of smoke plumes were observed largely between June and November, however smoke plumes 150 were continuously detected from Qayyara oil fields for a total of 225 days ranging from 13 th of June 2016 to the 27 th of March 2017. As a result of these event an estimated 1.33 million barrels of oil were burned (Bulmer, 2018;Tichý and Eichler, 2018). Residents south of the Qayyara oil fields were exposed for a number of 103 days to smoke plumes. Shortterm health effects were reported especially for patients with pre-existing respiratory conditions (Bulmer, 2018).
The Gulf of Sidraa has seen extensive episodes of smoke plumes as oil terminals at As Sidr and Ra's Lanuf, Libya 155 have been repeatedly set on fire over the course of a decade. These events have been capture by MODIS sensors through the last decade all the way since 2008. All events were characterized by dark plumes suggesting high contents of BC. On august 19, 2008 a tank fire erupted in Fiba tank farm at Ra's Lanuf after workers failed a maintenance operation (Piafom, 2018;The Telegraph, 2011). The fire lasted 9 days in which smoke plumes could be detected from the 19 th to the 22 nd . In March 2011 the terminal was struck by air artillery in the battle of Ra's Lanuf as the country was engaged in civil war (BBC, 2011;The 160 Christian Science Monitor, 2011;The Guardian, 2011). Smoke plumes were visible on the 12 th and 14 th . In December 2014 the tank farm at As Sidr oil terminal was struck by rockets as rebel fought to seize control of the city port. Seven storage tanks where engulfed in flames burning 1.8 million barrels of crude oil (Reuters, 2014;BBC, 2014;Al Jazeera, 2014).
Smoke plumes covered large areas of the Gulf between the 26 th and 30 th and could be seen as far east as Timimi and Crete.
In January 2016 both tank farms at As Sidr and Ra's Lanuf were struck by Islamic State militants. On January the 5 th As Sidr 165 https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License. suffered five tank fires while two tanks were hit at Ra's Lanuf amounting to 850,000 barrels of oil. A week later IS militants struck oil infrastructure connecting Ra's Lanuf terminal to other installations in the area (Tichý & Eichler, 2018;Tichý, 2019). A second attack at Ra's Lanuf tank farm was conducted on the 21 st (Business Insider, 2016). Throughout the month, extremely dense smoke plumes could be seen in the region from the 5 th all the way to the 23 rd . The most recent incident at Ra's Lanuf took place in June 2018 when rival armed groups clashed. The fire which started on the 14 th was contained at two 170 storage tanks before it was extinguished several days later (Reuters, 2018 Herein we will summarize the aerosol retrieval algorithms dark target (DT) over land and ocean Tanré et al., 1997;Levy et al., 2007aLevy et al., , b, 2009Levy et al., , 2013Remer et al., 2005Remer et al., , 2008Remer et al., , 2013, and Deep Blue (DB) over land (Hsu et al., 2004(Hsu et al., , 2006 with some emphasis on the atmospheric parameters used in the construct of lookup tables (LUT) and aerosol model selection as these properties/assumptions are crucial for proper AOD retrieval in oil smoke events.
The DT land algorithm is used over dark vegetated surfaces with low surface reflectance. DT makes use of the "VIS 190 to 2.1" relationship to distinguish surface contributions to the top-of-atmosphere (TOA) reflectance, as aerosols have a low absorbing and scattering effect in shortwave infrared (2.12 μm) band compared to the visible blue (0.47 μm) and red (0.66 μm) bands. After applying several pixel screening and selection techniques the algorithm chooses specific aerosol models and types based on seasonal and geographical characteristics. From here it determines aerosol related atmospheric parameters as a function TOA, surface and gas contributions to the apparent reflectance. To achieve AOD inversion these 195 parameters are matched to values within predetermined look up tables (LUT) as an attempt to describe the most likely aerosol conditions. The land algorithm uses five aerosol types composed of two or more models (fine or coarse), each with its specific aerosol optical properties. Primary products retrieved are AOD at 0.55 μm, fine model aerosol type (FMF) and spectral fitting error (ε).
The DT ocean algorithm works in much the same way as the land algorithm although it requires masking sediments 200 and filtering out strong glint areas. It uses the spectral dependencies of six bands, 0.55, 0.65, 0.87, 1.24, 1.63 and 2.12 μm for retrieving surface reflectance while it finds the exact match between the observed and the predetermined LUT reflectance values for the 0.87 μm band. It then attempts the best fit for the remaining bands. The 0.87 μm band is a good indicator of aerosol loading as it is less impacted by water radiance. Preliminary data such as AOD 0.55 μm, reflectance weighting parameter (η) at 0.55 μm, and aerosol effective radius (R eff ) are determined before the LUT inversion. DT ocean retrieves the 205 https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License. same products as DT land however it assumes aerosol properties based on a combination of two models, one fine model (four possible modes) and one coarse model (five possible modes). These models are combined by the weighting factor (η) such that the solution yields the lowest fitting error.
Aerosol properties within the LUT such as size distribution parameters, refractive indexes and SSA are crucial for proper aerosol typing and subsequent AOD retrieval. LUT information and model selection is critical for deriving other 210 optical properties, such as Ångström Exponent (AE) which can also be used to describe size distribution.
The DB algorithm was developed to retrieve AOD over arid, semi-arid and urban areas where surface reflectance values are higher than those over dark target regions. The principle behind the algorithm suggest that surface reflectance in these areas show higher values in red and near infrared bands and lower values in the blue band. The algorithm uses reflectance values form 9 bands through each step of the retrieval. After screening and pixel selection surface reflectance 215 values are determined based on three bands (0.412, 0.490 and 0.67 μm) using either a database or a dynamic approach. The approach is selected based on the normalized difference vegetation index (NDVI) and may be a function of region, season, scattering angle and land type. DB matches observed to LUT radiance values using a maximum likelihood method to determine the mixing ratio for a dust and a smoke model. The method retrieves two types of aerosol products: AOD and SSA from the dust model; AOD and AE from the smoke model. 220

CALIPSO aerosol data
The Cloud Aerosol Lidar with Orthogonal Polarization (CALIOP) onboard the CALIPSO satellite has been observing vertically distributed aerosol and cloud properties since 2006. CALIOP is an elastic backscatter lidar operating at two wavelengths (532 and 1064 nm) equipped with a polarization channel at 532 nm ).
Calibration is achieved through a molecular normalization technique for night-time measurements at 532 nm which is 225 subsequently the basis for daytime calibrations at both channels Kar et al., 2018). The latest CALIOP Version 4 data is significantly improved thanks to the refined calibration algorithms Kar et al., 2018;Vaughan et al., 2019).
CALIOP data requires several processing sequences, handled by different algorithms, to achieve the desired aerosol and cloud properties ). The first algorithm (selective iterative boundary locator -SIBYL) starts by 230 analysing calibrated level 1 data averaged horizontally (at resolutions from 0.33 km to 80 km) through the use of an adaptive threshold scheme establishing layer boundary limits . The next steps require the use of scene classification algorithms (SCA). Firstly, the cloud-aerosol discrimination (CAD) algorithm uses multidimensional probability density functions (PDFs) to distinguish clouds from aerosol layers (Liu et al., 2005(Liu et al., , 2009). The primary inputs from the CAD algorithm is later used for subtyping aerosol species throughout the troposphere and stratosphere (Omar et al., 235 2009;Kim et al., 2018). Finally, SCA uses the attenuated backscatter and volume depolarization ratios (both layerintegrated) to distinguish between water and ice clouds (Hu et al., 2009;Avery et al., 2020). To extract aerosol properties https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License.
(particulate backscatter and extinction coefficients, optical depth) the classified layer data is fed through several hybrid extinction retrieval algorithms (HERA) (Young and Vaughan, 2009;Young et al., 2013Young et al., , 2018. Lidar ratios are essential for calculating extinction coefficients, and throughout these sequences of algorithms lidar 240 ratios are selected in one of two ways. For unconstrained retrievals, lidar ratios are selected based on the aerosol subtype classification which is a function of surface type, location, particulate depolarization ratio, integrated attenuated backscatter (Omar et al., 2009;Kim et al., 2018). For each aerosol subtype a lidar ratio is assigned based on AERONET data, direct measurements and theoretical scattering calculations (Omar et al., 2009;Tackett et al., 2018). The second approach, known as constrained retrievals, is based on measured layer two-way transmittance (Young and Vaughan, 2009;Young et al., 245 2018). Selection between these two approaches is done based on scene complexity and feature classification (Young and Vaughan, 2009;Young et al., 2018). In most cases aerosol lidar ratios are determined using unconstrained retrievals (e.g., layers in contact with the surface) however constrained solutions are possible in certain situations (Young and Vaughan, 2009;Young et al., 2018).

Synergic approach 250
Figure 1 summarizes each step of the analysis starting from a collection of events identified in scientific literature and local news outlets. We chose a 12-year time frame in which both MODIS and CALIPSO retrieved a substantial amount of aerosol properties originating from oil smoke plumes. By using MODIS red, green and blue RGB images we can visually identify each plume (Table 1). The two MODIS sensors onboard Aqua and Terra were used as they possess a number of advantageous characteristics for plume identification and analysis: daily global coverage, good pixel resolution, algorithm 255 maturity, two retrieval windows, large data record (20 plus years of mission). CALIPSO data was used to compare plume AOD and to fill the gaps in MODIS data such as: plume thickness and elevation, scene classification, aerosol typing (lidar ratio and particulate depolarization ratio). Events identified in literature were selected for analysis based on plume dimensions and retrieval conditions. We selected plumes larger than 500 km 2 for statistical relevance as smaller plumes resulted in a low pixel count. Events were discarded if the atmospheric scene was predominantly cloudy, over 50% cloud 260 coverage. We deemed retrievals successful if they produce AOD pixel values, with some degree of variation (For at least 50% of pixels, AOD should vary resulting in value differences of at least 0.01), overlapping the plume area. We deemed unsuccessful retrieval if no AOD data was retrieved over the plume area (after cloud screening), or in cases where values of AOD were lower than 0.1 and the retrieved pixels were homogeneous (ex: over 90% of plume pixels with a fixed AOD value of 0.09 as seen in figure 6d). We used both successful and unsuccessful retrievals to highlight the capabilities and 265 limitation of the MODIS sensor. We used both MODIS Aqua and MODIS Terra AOD products, collection 6.1 (MODIS Atmosphere Science Team, 2017a, b), for successful retrievals. If products from both sensors were available for unsuccessful retrievals we selected only the most relevant scene. The algorithm for the AOD products was selected based on surface type (DT over ocean and land) and locations (DB over desert and arid areas) while aerosol properties were analysed only for successful retrievals. Based on our preliminary findings the DT ocean products were selected for successful retrievals while 270 https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License.
DT land and DB land products were used for unsuccessful retrievals. We took advantage of the higher resolution 3 x 3 km 2 level 2 AOD products for statistical relevance in successful retrievals over ocean. For unsuccessful retrievals we used the 10 x 10 km 2 level 2 AOD products (DB over desert and arid areas) and the 3 x 3 km 2 level 2 AOD products over land. The following aerosol properties were used in our analysis for successful retrievals: AOD at 0.55 µm, AE and R eff . To highlight the MODIS AOD product limitations we used AOD at 0.55 µm in scenes with unsuccessful retrievals. For successful 275 retrievals, we developed a plume averaging technique based on plume and background AOD values. Since both the RGB and AOD images show visible distinctions over the plume area, as seen in figures 2a and 2b., we constructed a plume edge based on AOD pixel gradient. The pixels selected to represent the plume edge share a unique AOD value which is different form the neighbouring background pixels by a value of at least 0.03. We averaged the AOD values from the plume area, within the plume edge, to obtain "plume average AOD". To achieve "plume specific AOD" values we subtracted a local averaged 280 background AOD value from the plume average AOD value. The result should represent the average AOD contribution of the oil smoke aerosols to the total column AOD, as seen in figure 2.c. Local background pixel count was chosen between 3 to 10 times the number of plume pixel count based on the local geography and meteorological conditions. A larger background pixel count is expected in scenes with low cloud coverage; over open water as seen in figures 3a and 3b. A smaller pixel count is expected in scenes with high cloud coverage over smaller gulf regions as seen in figures 3d and 3e. 285 Results and discussions for both successful and unsuccessful MODIS retrievals are presented in section 3.1.
One advantage of using CALIPSO data is that retrievals are done within 2 minutes following MODIS Aqua, for daytime retrievals, and thus plume and atmospheric conditions are similar for both sensors. Events were selected based on the plume cross section dimension. We used total attenuated backscatter (532 nm) values at single shot resolution to determine the extent of the plume cross section. Similar to MODIS data, we could identify a plume signature as backscatter 290 values in the cross section were larger than background levels. We chose plume cross sections larger than 5 km as level 2 data is averaged starting from this resolution. We used level 1B profile data (532 nm), standard version 4.10, both day and night time retrievals (Winker, David, 2016), to identify the plume signature. For daytime retrievals we used the Aqua MODIS RGB image prior to the CALIPSO retrieval for visual confirmation of the smoke plumes. For nigh time CALIPSO retrievals we used two consecutive daytime MODIS RGB images for an indirect "visual confirmation" of the smoke plume. 295 In these cases we used a MODIS RGB image the day prior to the CALIPSO retrieval and the first MODIS RGB image directly after the CALIPSO retrieval to assess the plume spatial continuity. We used Level 2 data (532 and 1064 nm), 5 km Aerosol Profile, standard version 4.20, both day and night time retrievals (Winker, David, 2018b) for statistical analysis of aerosol optical properties and comparison to MODIS AOD products. For aerosol typing and feature classification we used level 2 data of the Vertical Feature Mask product, standard version 4.20 (Winker, David, 2018c). We used the 5 km Aerosol 300 Layer product, standard version 4.20, both day and night time retrievals (Winker, David, 2018a) to assess the lidar ratio (532 and 1064 nm) if the plume was identified as a distinct aerosol feature. Focusing on the plume bins we used the extinction coefficient to backscatter ratio to determine the plume averaged lidar ratios if the plume was not identified as a distinct aerosol feature. We used a quality filtering procedure similar to the one described by Tackett et al., 2018. For cloud-free data we 305 used only Aerosol Profiles with a cloud-aerosol-discrimination (CAD) score of -100 < CAD score < -20. We discarded aerosol profiles in all range bins directly below any type of clouds as these profiles may be affected. Smoke plume profiles were discarded if the average plume height exceeded 4 km (above mean surface level) to prevent aerosols in contact with ice clouds and subsequent misclassifications and cirrus contamination. For the extinction filtering procedure, QC flag values not equal to 0, 1, 16, or 18 were discarded as low-confidence retrievals. Extinction samples where extinction uncertainty is 99.9 310 km -1 and retrievals in bins directly below these samples were rejected as these values may also be affected.
We selected the following level 2 aerosol data in the analysis: Particulate backscatter coefficients, extinction coefficient (532 and 1064 nm) and particulate depolarization ratio (532 nm). Based on the particulate total backscatter coefficient (532 nm) we defined the plume cross section as in each range bin, the plume values were at least 2 times higher than background values. Based on this parameter we also defined the plume thickness and elevation. We defined the 315 background region as the same number of bins and of the same height as plume bins, situated either upwards or downwards of the plume location (depending on plume trajectory), following the overpass trajectory. The extinction values (532 and 1064) in plume bins were used to integrate the AOD for each plume profile. The resulting AOD values were averaged to obtain a plume mean AOD value. For both plume cross section and background area we averaged the aerosol optical properties (particulate backscatter, extinction and depolarization ratios). To assess the possibility of aerosol typing we 320 averaged the plume lidar ratios and calculated the Ångström exponent using the plume mean AOD values (532/1064 nm).
Averaged data was compared between plume cross sections and local background values to assess CALIPSO feature detection and aerosol typing algorithms. We compared our data between events to determine how lidar ratio solutions affect aerosol properties estimates. Based on particulate depolarization ratios and lidar ratios we assessed the possibility of oil smoke aerosol typing either as a distinct class or as part of the predetermined CALIPSO aerosol types. 325 Since both MODIS and CALIPSO identified plume specific AOD and AE values we compared our results to determine each methods strengths and weaknesses. We also identified one AERONET event which was also used to compare plume values. We compared our results to other studies done on oil smoke plumes in which the same aerosol properties were determined by means of ground based or flight measurements. Our conclusions reflect the nature of these smoke plumes and the implications they have on current satellite retrieval capabilities.

340
We selected a successful retrieval to better describe the method used for our analysis. Figure 2 shows the case at Ra's Lanuf and As Sidr tanks farms which caught fire on the 5 th of January 2016 and burned throughout the 6 th and 7 th . The retrieval in the images was taken on the 6 th by MODIS Aqua at 12:05 UTC. mixing ratio with the local background aerosols. In this study, we focused our attention on the plume areas where heavy concentrations of aerosol are present while discarding retrievals done at the edges of the plume where background aerosol may have a large influence on the values retrieved. Thus figure 2b was constructed based on the AOD (0.55 µm) retrieval and a plume edge selection technique. To determine the plume edge we constructed isolines of AOD values from the 350 retrieval pixels. The 3×3 km product is better suited for determining the AOD gradients and thus was selected over the standard 10×10 km product. Figure 2b shows higher values of AOD in the selected plume area as opposed to the local background levels. At the time of the retrieval, we can observe two distinct plumes of smoke, a thin plume originating from the As Sidr and the main plume (within the black contour, Fig. 2b) originating from Ra's Lanuf. Since the As Sidr plume did not meet the selection criteria the analysis is made for the main Ra's Lanuf plume. To further discriminate between plume 355 and background AOD values we averaged all non-plume pixel values, over water, within the gulf region without considering the pixels of the plume. Then, the averaged AOD values were subtracted from each pixel of the plume AOD values to determine the overall plume contribution. Consequently, Figure 2c illustrates the plume specific AOD gradient.   No.10 was operating 24 oil wells and 4 gas wells of which all 4 gas wells were involved in the fire. The resulting smoke plume can be seen from MODIS true colour images exhibiting a light grey, whitish colour. The plume albedo is more similar to biomass smoke and less similar to black smoke plumes as seen from the events presented in this study. Thus, a correct LUT construction and matching of the observed atmospheric properties should take into account the plume's SSA. In figure   3a we observe the plume's specific AOD and AOD gradient computed from MODIS Aqua retrieval. As seen in Tables 2 and  380 3, mean AOD plume specific values fall between two and three times higher than local background values. Within a threehour window between the two retrievals (Terra vs Aqua) we can observe a slight drop in AOD values as the plume is being dispersed, however these small changes may also be attributed to other factors such as the satellite sensor calibration and retrieval geometry. As opposed to Platform No.10 fire, the Deepwater Horizon fire (Gulf of Mexico) was mainly fuelled by oil. This is evident in the plume albedo from MODIS true colour images. The gradient in figure  Nippon fire produced AOD values two times higher than the background values, within the first hours of the event, after which values declined slowly at the time of Aqua overpass. Figure 3c shows a plume specific AOD values ranging from 0.06 390 to 0.23. Figure 3c shows AOD values as high as 0.24 over the average AOD background level for the plume originating at Ra's Lanuf. During this event, an increase in AOD levels between the two retrievals was observed as the fire spread to several oil tanks. The AOD gradient, in figure 3c, shows the largest values at the centre of the plume where aerosol mixing is expected to be lower. Two weeks later the same tank farm injected thick plumes, figure 3d. over the Surt district in Libya.
The plume section over the Gulf recorded AOD values twice as high as the background level however the net contribution 395 amounted, on average, to a value of 0.10. Plumes from As Sidr are visible in Figs. 3e, 3f, and 3g. This event was captured in multiple days while the fire engulfed several oil tanks and subsequently injecting higher amounts of aerosols in the region.
Depending on the local background levels, average plume specific values ranged from -0.03 to 0.15. Negative values can be explained by the presence of dust and marine aerosols in the atmospheric background. This is especially evident in figure 3g where high background levels were registered in the Gulf of Sidra while lower levels were seen off the shores of At Tamimi, 400 600 km NE of As Sidr. The Puerto Sandino fire produced a thinner layer of smoke which was evident in the AOD retrieval with values close to the background level. As seen in figure. 3h the smoke plume is largely dispersed and thus exhibiting lower levels of AOD as opposed to 3 hours earlier (see table 2  We computed the Ångström exponent (0.55/0.86 µm) for plume and background values as previously described. 415 The outlier case from SOCAR's Platform No. 10 fire shows a slight difference in AE plume values with respect to the background aerosols. As seen in table 2 and table 3 table 2 and table 3. While AE plume values are generally low, these extremely low values may not be primarily a direct result of particle size distribution. MODIS uses spectral reflectance relations to determine AOD and subsequently AE levels. While other types of aerosols have a varying spectral reflectance signature, heavy concentrated black carbon exhibit a flat and linear signature that result in low spectral reflectance values (Johnson et al., 1991;King, 1992;Pilewskie and 430 Valero, 1992;Soulen et al., 2000).  After analyzing plume values seen in figures 2 through 5 and tables 2 and 3 we can determine that MODIS algorithm over ocean is able to detect large changes in AE and R eff values. Aerosol load however is only slightly impacted by the presence of heavy smoke and thus to validate these results, other datasets such as CALIPSO and AERONET data were 465 used.   This may be a result of the algorithm not finding a proper fit to the LUT. In these cases, radiance, surface and TOA reflectance values are lower as a result of the smoke plumes low transmission and albedo. Thus, predetermined values within the LUT may not have included such "extreme" cases. Figure 6d shows the event at As Sidr also captured by the ocean algorithm and analysed in the previous section (successful retrievals). In this case we can distinguish the plume from the 480 RGB image over the Gulf of Sidraa while also observing AOD values over land where the smoke plume drifted E-NE towards the island of Crete. However, there seems to be no distinguishable AOD gradient within the plume. A further inspection suggested that all pixels showed values of 0.095 which suggest that the lower radiance values did not match well with pre-existing LUT values. Consequently, the region is classified as "clean atmosphere" and thus, a unique AOD value is assigned to all the pixels. Conversely, the ocean algorithm retrieved AOD that varied between 0.1 and 0.37. The same 485 situation emerged from the Jaipur event (figure 6e) where plume pixels did not present any AOD gradient with only pixel AOD values of 0.095. Background values are much higher averaging 0.3 while biomass burning in the area recorded values of AOD over 1. It is safe to say that AOD values for black smoke aerosol are not correctly retrieved over land surfaces due to lower atmospheric transmission all values are either incorrectly selected by LUT or are assigned to a "clean atmosphere" value. Figure 6f and 6g show the unsuccessful retrievals for San Martín Texmelucan de Labastida of Puebla, Mexico and 490 Vasylkiv, Kyiv Oblast, Ukraine. It was expected for both cases that larger values of AOD should be observed over plume pixels with lowest albedo. However, in these regions the land algorithm produces fill values and thus, no AOD value could be retrieved. In figure 6f we can identify higher AOD values west of the black plume area that correspond to local biomass burning events and are successfully retrieved since these atmospheric conditions are included in the LUT. It is also obvious that at the edge of the black plume, the algorithm retrieved a negative value which leads us to believe that values lower than -495 0.10 were identified within the plume. This may only be the result of unmatched MODIS observed to LUT values. Since these heavy smoke plumes are the result of extreme scenarios they are rarely observed and may not end up being a subject of research. Thus, we believe there are no cases within the LUT values describing extremely low atmospheric transmission and radiance values, highly absorbent aerosol, low SSA and low reflectance values over a large spectral range including MODIS bands 1 through 7. 500

CALIPSO Data
The event at Ra's Lanuf and As Sidr, 6 January 2016, was also captured by CALIPSO lidar measurements as CALIPSO overpass matched a cross section of the plume area. Figure 7a shows this overlap in near-real time as CALIPSO succeeds Aqua within a 2 minute time frame. Within the 15 km plume cross section we selected a particulate backscatter 505 coefficient profile for reference, figure 7b, and based on this parameter we determine plume elevation and thickness. The average plume thickness was approximately 920 m ranging from 2700 m to 3300 m. The entire plume cross section is presented in figure 8. We observe the main plume from Ra's Lanuf elevated between 2400 m and 4200 m. values at 1064 nm measured 0.17 km -1 sr -1 . Average extinction coefficient values at 532 nm were measured at 1.65 km -1 while the 1064 nm channel yielded a value of 1.55 km -1 . This event is an example of an opaque aerosol layer, were the lidar did not penetrate the plume up to the sea surface over the Gulf of Sidra. Table 6 shows a lidar ratio of 109 sr at 532 nm and 86 sr at 1064 nm. These values are larger than the CALIPSO V4 aerosol subtype values for: elevated smoke 70 ± 16 (532 nm) and 30 ± 18 (1064 nm); Polluted continental/smoke 70 ± 25 (532 nm) and 30 ± 18 (1064 nm) (Kim et al., 2018). The 515 particulate depolarization ratio for the Ra's Lanuf plume was 0.11 that corresponds to moderately depolarizing smoke (Kim et al., 2018). Figure 8c shows the CALIPSO feature classification while figure 8b shows the aerosol typing results.  Judging from these images and from the average CAD score of -48, the smoke plume represented a mixed feature of cloud and aerosols. This is to be expected as water vapour and particulate matter are primary components in emissions of petrochemical burnings (Johnson et al., 1991;Ferek et al., 1992;Daum et al., 1993). Cloud formations on top of oil smoke plumes have been observed in other instances (Johnson et al., 1991), as seen in figure 9, a phenomenon which hinders AOD retrievals in both passive and active sensors.  The current version of the Vertical Feature Mask gives a mixed result to aerosol typing comprised of dust, polluted dust and smoke aerosols for this oil smoke plume. Table 5 shows large values for AOD ranging from 1.52 (532 nm) to 1.43 (1064 nm). The low background AOD values indicate a clean atmospheric scene, thus the large AOD values are considered to be related to the smoke plume. We also computed the AE (532/1064) which resulted in 0.09, indicating the presence of 540 coarse particles. It should be mentioned that this event was an optimal case for constrained lidar ratio retrieval since the feature was surrounded by clear air. However, due to the fact that the plume feature was completely opaque, the lidar ratio could not be obtained from a constrained solution via the two-way layer transmittance. Within the 12-year period we identified two more cases in the Gulf of Sidra. The event at As Sidra on 29 December 2014 is an example of an unconstrained lidar ratio solution. In this case, as with all the remaining cases, CALIPSO nighttime overpasses do not allow for a direct comparison with MODIS data. The plume at As Sidr can be observed in figure 10a.
Off the coast of As Sidr a distinct feature above the sea surface reaching 650 m in altitude is observed. Figure 10b show a 555 particulate backscatter profile where we can distinguish a plume thickness of approximately 240 m. This event was smaller in magnitude with respect to the Ra's Lanuf event where multiple storage tank fires contributed to the same plume mass. In this case, CALIPSO overflew much closer to the tank farm also resulting in a narrower plume cross section. The SIBYL algorithm failed to detect the plume area as a specific feature and thus, level 2 products were averaged over a larger 20 km area, as opposed to the 5 km averaging resolution, in which background aerosol levels were significantly lower. This was 560 expected as the plume cross section measured approximately 3 km, as seen in figure 10c and d. Consequently, CALIPSO identified dusty marine aerosols within the Gulf region as evident from the lidar ratio of 37 sr (both 532 and 1064 nm) within the atmospheric layer (Kim et al., 2018). The plume particulate depolarization ratio value, 0.12 was close to local background values of 0.18 and similar to the values from Ra's Lanuf event, 0.11. Although this value is similar to the previous event, it should be mentioned that this value is indicative of the larger 20 km feature and not specifically to the 565 https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License. smaller 3 km plume feature which was not successfully identified by the SIBYL algorithm. The local scene over the Gulf region was fairly clean judging by the low AOD values seen in table 5. The low plume AOD value can be a result of the larger averaging resolution within the background layer. Another reason for low AOD retrievals was proposed by Jethva et al. (2014) (Jethva et al., 2014) and was shown in Deaconu et al. (2017) (Deaconu et al., 2017). In case of optically thick aerosol layer, the sensitivity of the backscattered signal would be reduced or lost because of the strongly attenuated two-way 570 transmission. As a result, the operational algorithm may position the base of the aerosol layer higher in altitude, thus underestimating the geometrical thickness of the aerosol layer and consequently the AOD. The selection of an inappropriate aerosol lidar ratio might also contribute to the underestimation of the AOD. AE values were also relatively low, 0.12 however these are again not indicative of plume optical properties as they are for the larger background layer. In any case, computing AE based on low AOD values may not be a good estimate for local aerosol particle size. Figure 10d shows the 575 same plume composition as in the case of the Ra's Lanuf plume. The vertical feature mask shows a mixt feature of clouds, aerosols and low confidence aerosol. Figure 10c, shows the aerosol classification within the aerosol layer. Judging by these results the aerosol layer reaching up to 1.5 km above mean sea level was classified as a mixture of dust, polluted dust and smoke. The plume cross section was successfully identified as smoke aerosols however the background layer situated north of the plume was not typed as dusty marine, as seen from the lidar ratio of 37 sr.  Another event in the region was captured inland as the plume cross section was identified 170 km south of the Ra's Lanuf tank depot. Figure 10e shows a particulate backscatter profile through the plume center describing a fairly inhomogeneous mass of smoke particles. The main plume was concentrated between 500 and 900 meters however lower concentration may have been mixed with local dust particles all the way up to 1500 meters. Figure 10f shows the extent of 590 the plume as it travelled southwards inland. As was the case of the previous event, plume lidar ratios were determined by an average plume particulate backscatter (532 nm) measured 0.007 km -1 sr -1 while the 1064 nm channel measured 0.008 km -1 sr -1 which showed an increase of 6 to 9 times larger than the local background values. This large difference is also evident from the AOD values since averages of 0.16 (532 nm) and 0.15 (1064 nm) were 5 to 9 times higher than the local background levels, however these plume values were directly influenced by the lidar ratio solution. Thus, a constrained solution may have resulted in larger values since smoke LR values are generally higher than polluted dust values (Kim et al., 2018). AE 600 values remain low 0.05 and similar to the previous event at Ra's Lanuf suggesting a coarse dominant aerosol mixture, while background values were more indicative of polluted dust, averaging 0.71. Unlike the As Sidr event, this plume was correctly identified by the SYBIL algorithm thus resulting in plume AOD values higher than background values. Similar to the previous events the feature classification algorithm shows a mixture of clouds and aerosols in the plume bins. This composition is evident in figure 10h and may affect the retrieval of smoke optical properties. Figure 10g shows the local 605 aerosol layer up to approximately 2 km above local ground level. A significant portion of plume bins were discarded as cloud features even though the CAD score of -99 indicated high confidence aerosols.

615
The events at the Qayyara oil fields in Norther Iraq were captured by CALIPSO in three distinct cases. In all three cases CALIPSO overpassed within less than 35 km southwest from the well fires. The first case was captured on 1 July 2016 as seen in Figure 11a. The plume can be observed just above the local surface travelling southwards and spanning over 100 km. Figure 11b shows the plume elevation which was identified from approximately 250 m (local surface elevation) to around 400 m above sea level averaging 150 m in thickness. The plume particulate backscatter coefficient values ranged 620 from 3 (532 nm) to 5 (1064 nm) times higher than local background values. Since the plume was identified within the PBL, a https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License. lidar ratio of 44 sr for both channels was assigned indicating the presence of dust aerosols. Particulate depolarization ratio values were higher than the previous cases, on average 0.27 for plume and 0.25 for background, suggesting dust dominance.
AOD values observed in plume the bins are 2.7 times higher than the local background. As with the previous cases the AE values remain consistently low 0.39 indicating the presence of coarse dominant aerosol mixture. Figure 11c shows the 625 aerosol classes identified in the local scene. The results indicate a dust dominant mixture with no obvious smoke signatures.
As opposed to the previous cases in Norther Libya, the vertical feature mask did not indicate the presence of clouds features in the plume bins. This may be caused by low relative humidity levels and high percentage of dust particles within the mixt feature.
The second case was observed on 17 July 2016 (figure 11d) locating the plume bins closer to the burning oil wells. 630 The plume was much narrower than the previous case however particulate backscatter coefficient values remained similar for both cases, 0.007 km -1 sr -1 (532 and 1064 nm). Figure 11d shows the plume top height reaching 450 m above sea level with elevated particulate backscatter levels in the first 3 bins above local ground level, figure 11.e. As with the previous case the plume was contained within the PBL and a lidar ratio of 44 sr for both channels was assigned indicative of dust particles.
The particulate depolarization ratio showed large values, i.e., 0.32, which indicate a dust dominated scene. Plume AOD 635 values were on average 2 to 3 times higher than local background values while AE remained low at -0.03 similar to the previous events. The feature classification algorithm did not identify a cloud presence in the plume bins, with figure 11.f showing the presence of dust particles.
The last event retrieved was very different to the previous two events at Qayyara. On 21 October, 2016 the black smoke plumes from the oil fields were mixed with a thick SO 2 plume as a result of the Islamic State setting fire to the Al-640 Mishraq Sulphur plant situated NNE of the burning oil fields. Figure 11g shows the oil smoke plume over the Qayyara oil fields with average particulate backscatter coefficient (532 and 1064 nm) values of 0.14 km -1 sr -1 . Judging by figure 11h and 10i we can deduce that a mixed aerosol layer was distributed between 1500 and 2600 m. This mixed layer was also suggested by (Kahn et al., 2019) who used MISR Active Aerosol Plume-Height (AAP) to establish the SO 2 and oil smoke plume elevation on the same day (Khan and Zhaoying, 2020). Regarding the plume close to the surface, which we attribute 645 to be oil smoke, the particulate backscatter (532 nm) plume values were 3.9 times higher than local background values. The smaller oil plume was observed, on average, 120 m above the local surface levels. Lidar ratios of 47 sr and 45 sr (532 and 1064 nm) for the PBL containing the plume suggested the presence of dust particles. The plume particulate depolarization ratio however revealed lower values of 0.15 versus 0.22 for the background layer suggesting the mixing with smoke aerosols. In this case the smoke plume increased the overall background values of AOD by a factor of 4. Since the 650 unconstrained solution for the lidar ratio did not indicate the presence of the smoke plume, it is safe to say that a constrained, smoke-like solution would have resulted in higher AOD values. As with all cases previously described a low AE value of 0.15 indicated the presence of coarse dominant dust and smoke particles. The aerosol subtyping algorithm did not distinguish between oil smoke, SO 2 and dust particles revealing no obvious smoke signatures. This is evident in figure 10.i as the local scene is seen dominated by dust and polluted dust. Based on CALIPSO measurements, the smoke backscatter and extinction coefficient ranged from 2 to 9 times higher than background levels. In four out of six cases, particulate depolarization ratio revealed values between 0.11 and 0.15 resembling moderately depolarizing smoke while larger values in two cases were mostly due to the presence of dust particles in the local atmospheric scene. Apart from one case, all lidar ratios were assigned to a constant value as the plume resided in the PBL. The opaque feature measured high lidar ratios of 109 sr (532 nm) and 86 sr (1064 nm) that resembles 660 smoke lidar ratio found in literature (Giannakaki et al., 2016;Haarig et al., 2018). of events it is difficult to assign a separate aerosol type for these oil smoke plumes. However valuable information regarding size distributions, particulate depolarization ratio and to some extent lidar ratio can be retained from this study. It should be mentioned that these values reflect smoke plumes located very close to the fire sources and thus present low mixing ratios with other local aerosols.

AERONET case study 675
As discussed in the introduction section, oil smoke plumes have been rarely observed using ground based remote sensing instruments such as AERONET sun photometers. We used AERONET version 3 direct sun data to assess the presence of oil smoke plumes. Only one study was found in scientific literature (Mather et al., 2007) which measured aerosol properties of the Bouncefield plume at two distinct locations. Here we identified the smoke plume resulting from naphtha tank fires in Vasylkiv, Kyiv Oblast, Ukraine on 9 June 2015. The smoke plume was also captured in RGB images as seen in 680 figure 6g. Figure 12a shows the distinct signature of the oil smoke plume as AOD values increased significantly in all wavelengths. Figure 12c is a good indication of the increasing particle size with respect to the other days observed as well in MODIS and CALIPSO data. Figure 12d shows the daily evolution of AE with values between 0.45 and 0.9 for the time frame in which the plume was observed. Figure 12b shows AOD values rising as the plume was travelling NE over Kiev.
The AERONET station in Kiev is situated approximately 35 km NE of the Vasylkiv tank farm. The peak of the plume was 685 detected at 9:45 UTC when the AOD was 0.68 at 500 nm. Unfortunately, no inversion products coinciding with direct sun measurements were available as the Kiev sky was partially cloudy at the time.   Table 7. Oil smoke optical properties from ground based and flight measurements along with the scientific reference.

Aerosol properties
Reference Comments Gulf War smoke plumes AOD: 0.82 -1.92 (Pilewskie and Valero, 1992) Visible optical depths AOD: 1.5 ; AE: 0.7 (Nakajima et al., 1996) 0.5 µm AOD: 2.0 -3.0 (Hobbs and Radke, 1992) Visible optical depths PDR: 0.14 (Okada et al., 1992) 532 nm (lidar)  (Ross et al., 1996) 532 nm (lidar) Extinction coefficients: 0.06 -1.30/km; 0.06 -1.60/km AOD: 0.05 -1; 0.05 -1.20 (Laursen et al., 1992) 1064 and 532 nm (lidar) Extinction coefficients: 0.5 -10/km (Weiss and Hobbs, 1992) 538 nm (photometer) AOD: 0.5 (Johnson et al., 1991) 11 µm Bouncefield smoke plume Near source: AOD 0.4 -1.4; AE 0.42; R eff 0.45 -0.85 µm Distant location: AOD 0.3 -0.7 AE 0.09; R eff 0.88 -1.40 µm (Mather et al., 2007) 440, 675, 870, 936 and 1020 nm (sun photometer) The results presented in this study show a wide range of values which are attributed to a multitude of local factors 695 such as: background aerosols, burning rates, weather conditions, fuel type, time of retrieval, and local geography. Other factors can be attributed to the different types of methods and algorithms used to retrieve aerosol specific data. MODIS data showed relatively low values of plume specific AOD ranging from -0.04 to 0.16. The only event which was captured by both MODIS and CALIPSO retrievals showed a large level of discrepancy. In particular, the Ra's Lanuf event on 6 January 2016 showed average column AOD values of 0.21 over the plume area with a maximum pixel value of 0.32 (550 nm). In contrast, 700 CALIOP measurements revealed average plume AOD values of 1.52 (532 nm) which where plume specific as no other extinction values were detected beneath or above the plume through the troposphere and stratosphere in the local scene. In the remaining 5 cases, CALIPSO retrieved AOD values ranging from 0.02 to 0.16 for average plume thickness ranging from 0.120 to 0.375 km. While these values more closely resemble the successful MODIS retrievals one should restrain from a forward comparison. A reason is that neither of these events were analysed by both sensors since MODIS did not 705 successfully retrieved AOD values over land. This coupled with high levels of uncertainty surrounding MODIS LUT values and CALIPSO lidar solutions suggest the need for a more in depth analysis. The one case seen by an AERONET sun photometer indicates AOD values ranging from 0.3 to 0.68 (500 nm) however the satellite images suggest that these values were not indicative for the main plume which most likely did not reached Kiev. Nevertheless, MODIS did not successfully retrieve any AOD values from this event or any other over land while for other events over ocean it did not yield such high 710 AOD values. It is safe to say that MODIS AOD retrievals for oil smoke plumes may not produce satisfactory results since the predetermined LUT values may not contain similar events to the ones described in this study. The CALIPSO AOD measurements are directly influenced by the lidar ratio. For these specific events, a correct estimate of lidar ratio is very difficult to achieve based on unconstrained solution. On one hand these lidar ratios are not directly measured. On the other hand, lidar ratio for oil smoke plumes may exhibit a different behaviour, considering the high BC content, different from 715 biomass smoke or smoke/polluted continental aerosols. In cases of "clean" atmospheric conditions a constrained solution may result in better AOD estimates however these conditions are rarely achieved, with less than 0.01% of all aerosol layers detected . AE values seem to be more consistent between MODIS, CALIPSO and AERONET as all https://doi.org/10.5194/acp-2021-790 Preprint. Discussion started: 13 October 2021 c Author(s) 2021. CC BY 4.0 License.
techniques suggest the presence of coarse aerosol mixtures, however in conditions with low AOD values, one should restrain from direct comparisons. 720 Table 7 lists the oil smoke optical properties from different studies that on average indicate much larger AOD values. It should be mentioned that AOD values from the Gulf war smoke plumes are larger for the most part due to the magnitude of the event. These measurements describe super composite plumes resulted from o large number of well fires and pool fires. An event that more closely resembles the events in this study was analysed by Mather et al. (2007) who retrieved AOD values of 0.3 to 0.7 50 km away from the oil depot (Mather et al., 2007). These values are similar to the 725 AERONET values from Kiev presented in this study as in both cases AOD was retrieved using sun photometers and the Kiev AERONET station is located at approximately 20 km from the oil depot. However, the Bouncefield event was significantly larger than the Vasylkiv event. CALIPSO AOD measurements from Ra's Lanuf (January 6, 2016), are similar to other studies described in table 7 while AE values are in general in good agreement with our case studies. Judging by the results seen in figure 8.a, an aerosol feature exhibiting such heavy attenuation in the layers directly beneath would most 730 likely yield higher AOD values. R eff values from MODIS also show the presence of larger particles analogous to Mather et al. (2007). Particulate depolarization ratio for four out of six cases reflects the values shown by Okada et al. 1992 indicating that oil smoke particles are moderately depolarizing. The opaque feature indicated high lidar ratio values (109 sr) much larger than (Ross et al., 1996). This may be a result of plume dimensions as the Ra's Lanuf plume cross section was much larger and thicker than the plume described in Ross, 1996. One would expect large lidar ratio values to be the result of the 735 highly absorbent nature of these smoke plumes (high percentage of black carbon, high plume homogeneity and low mixing ratio) leading to larger extinction values. CALIOP extinction coefficient averaged for each individual plume in this study ranged from 0.10 to 1.65 km -1 (532 nm) and 0.10-1.55 km -1 (1064 nm) with a maximum bin value of 9.6 km -1 . These values agree well with (Weiss and Hobbs, 1992) and (Laursen et al., 1992).

Conclusions 740
In this study, we examined oil smoke plumes derived from 30 major industrial events within a 12 year period. To our knowledge this is the first study that utilized a synergetic approach based on satellite remote sensing techniques. The MODIS ocean algorithm was successful in retrieving aerosol properties which, on average, ranged from -0.06 to 0.16 for plume specific AOD, -0.18 to 1.25 for Angstrom exponent and 0.29 to 1.73 µm for effective radius. Apart from the SOCAR event, all the remaining smoke plumes exhibited AE values lower than 0.74 suggesting that the smoke plumes were coarse 745 mode dominant. CALIPSO measurements showed values of plume AOD ranging from 0.02 to 0.16 (532 nm) and 0.02 to 0.15 (1064 nm) except for one event where AOD values reached 1.52 (532 nm) and 1.43 (1064 nm). AE values ranged from -0.03 to 0.39 which agree with MODIS. A large discrepancy can be found in one event (at Ra's Lanuf, Lybia, on the 6 th of January 2016) where CALIPSO measured column integrated AOD values 5 times higher than MODIS. For this event CALIPSO produced a retrieval that resulted in high lidar ratios of 109 sr (532 nm) and 86 (1064 nm). The high concentration of water vapour emitted by the oil fire may have contributed to instances of small cloud formation above the smoke plume and thus contaminating the retrievals. Typically, lidar ratio ranged from 37 to 55 sr (532 nm) and 37 to 48 sr (1064 nm) however these unconstrained solutions were indicative of the local aerosol scene and not directly measured. Particulate backscatter coefficient values ranged from 0.002 to 0.015 km -1 sr -1 (532 nm) and 0.002 to 0.017 km -1 sr -1 (1064 nm).
Particulate extinction coefficient values ranged from 0.10 to 1.65 km -1 (532 nm) and 0.10 to 0.155 km -1 (1064 nm). On 755 average backscatter and extinction coefficient values were 2 to 9 times higher than the local background. Particulate depolarization ratios ranged from 0.11 to 0.15 for four out of six cases while the remaining cases were 0.27 and 0.32. We suspect that this discrepancy in the two cases at Qayyara are a result of dust aerosols presented in the smoke plumes. The values presented agree with similar studies that used ground-based and airborne measurements. We believe that MODIS gives a conservative estimate of the plume AOD since MODIS algorithms rely on general aerosol models and various 760 atmospheric conditions within the look-up tables which do not reflect the highly light absorbent nature of these smoke plumes. Furthermore, the spectral reflectance relationship used by MODIS algorithms may hinder most retrieval attempts as thick black plumes exhibit a distinct spectral signature. CALIPSO measurements are heavily dependent on unconstrained solutions, which in turn do not reflect the oil smoke plumes. Thus, we also believe that the AOD values based on CALPSO measurements are conservative in nature since strong absorbing smoke would yield larger lidar ratios and AOD values. We 765 stress the need for further lidar measurements of oil smoke plumes since, based on this study we cannot conclude whether these aerosols belong to a different smoke subtype. Future space-borne lidar mission such as EarthCare (Illingworth et al., 2015) will provide direct measurements of lidar ratios and the possibility of better AOD estimations with regard to these types of events. Based on this study we concluded that the MODIS land algorithms are not yet suited for retrieving aerosol properties for these types of smoke plumes due to the highly light absorbing properties of these aerosols. This study has 770 shown a novel method of oil smoke plume identification and analysis which does not require, in some cases, perilous field work. We believe that these types of studies are a strong indication for the need of improved aerosol models and retrieval algorithms. For these types of aerosols, better AOD estimates are important for both air quality and climate change implications.

Code availability 775
Not applicable.