Retrieval of Land Surface Temperature from Earth Observation Satellites for Gas Flaring Sites in the Niger Delta, Nigeria

This research investigates the recording of Land Surface Temperature (LST) by Earth Observation (EO) Satellites for four gas flaring sites in Rivers State, Nigeria. Six Landsat 5 Thematic Mapper (TM) and Eleven Landsat 7 Enhanced Thematic Mapper Plus (ETM+) from 17 January 1986 to 08 March 2013 with < 5% cloud contamination were considered. All the sites are located within a single Landsat scene (Path 188, Row 057). Dark Object Subtraction (DOS) method and Atmospheric Correction Parameter (ATMCORR) Calculator were used to obtain atmospheric correction effects parameters for multispectral and thermal bands [Upwelling radiance (Lu), downwelling radiance (Ld) and transmittance (τ)] of Landsat data respectively. The emissivity (ε) for each site is estimated by using standard values for determined land surface cover from Look Up Table (LUT). The correction obtained from DOS method was applied to the computed reflectance to get the atmospherically corrected reflectance that was used for the classification of land cover. The Lu, Ld and τ obtained were applied to the calibrated at-sensor radiance band 6 (high gain) data to compute the surface-leaving radiance (Lλ) with the ε values obtained for each site. The Planck equation was inverted using the calibration constants to derive LST. Six range of LST values were retrieved for each flaring site, with Bonny Liquefied Natural Gas (LNG) Plant recorded the highest LST (345.0 K) and Umudioga Flow Station with the lowest (293.0 K). LST retrieved from both sensors for the flare hotspots are the highest values compared to other locations within the processing sites, which was clearly shown through Geospatial Information System (GIS) spatial analysis and the transects plots. Furthermore, the closer is the distance to the flare, the higher is the temperature and vice versa. Based on these results, it can be concluded that satellite based sensors, such as Landsat TM and ETM+, have the ability to record LST at gas flaring sites in the Niger Delta.


Introduction
Remote sensing of Land Surface Temperature (LST) is an important research area due to its diverse applications such as detection of gas flares or forest fires [1,2], land use and land cover change (LULC) [3], vegetation monitoring [4], climatic change analysis [5], geothermal area detection, weather prediction and analysis of energy and matter exchanges between the atmosphere and surfaces [6,7], hydrological circle [8]. LST variations in space and time, measured by satellite remote sensing, are used for the estimation of a multitude of geophysical variables, such as evapotranspiration, vegetation water stress, soil moisture, and thermal inertia [9][10][11].
Due to the increasing recognition of the importance of LST, methods for its estimation from space have continuously been developed [12]. For example, sensors, such as the Advanced Very High Resolution Radiometer (AVHRR) [13] and the Moderate-resolution Imaging Spectroradiometer (MODIS) [14], have provided public domain global thermal data twice daily, using two longwave infrared (LWIR) bands. Advanced Along Track Scanning Radiometer (AATSR) also provide thermal data with two LWIR at every 35 days [15]. Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) provide thermal data using just one LWIR band 6, while Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIR) provide thermal data with two LWIR bands 10 and 11 [16]. The three Landsat sensors has a higher spatial resolution but with a 16 day temporal resolution. Furthermore, [17] integrated Sentinel-3 Sea and Land Surface Temperature Radiometer (SLSTR) data and Visible and Infrared Scanner (VIIRS) data to record LST for investigating gas flaring activities; [18] used MODIS data for the retrieval of LST for the characterization of gas flaring from the space; and [19] compared MODIS data, Landsat 8 OLI/TIR data, VIIRS data and Spining Enhanced Visible and Infrared Imager (SEVIRI) data for retrieval of fire radiative power (FRP) for assessing industrial gas flaring output. Since satellite remote sensing provides a repetitive synoptic view in short intervals of the Earth's surface, it is a vital tool for monitoring LST [20,21]. However, there are some challenges that affect the accuracy of the derived LST values obtained from satellite data. For example, the upwelling thermal radiation measured by satellite systems which is used as a substitute for estimating the LST is affected by atmospheric constituents before reaching the sensors resulting in inaccurate LST estimates if the atmosphere is not correctly accounted for [22]. Several approaches have been developed to handle the problem of atmospheric correction for thermal bands. Among these are the single or multi-channel (split window) [23] and dual angle approaches [1]. A variety of split-window methods have been developed to retrieve LST from National Oceanic and Atmospheric Administration (NOAA)/AVHRR data [24]. The split-window LST method utilizes the differential absorption in adjacent thermal band to correct the atmospheric effects [25,26]. [27] propose a multi-band algorithm to retrieve LST and land surface emissivity from MODIS, which is only influenced by the surface optical properties and the ranges of atmospheric condition. The accuracy of these two algorithms is less than 1°C [28,29]. [30] also propose an algorithm to retrieve temperature and emissivity from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). The accuracy of the results dependon the empirical relationship between emissivity values and spectral contrast, compensation for reflected sky irradiance, and ASTER's precision, calibration, and atmospheric compensation. However, the accuracy of most algorithms is very high but they still need to make some assumptions regarding prior knowledge of atmosphere (especially water content). Owing to different considerations of the atmospheric effect on the radiation transfer through the air, the prior knowledge required is different [24]. [31] and [32] concluded that including column water vapour in the split-window algorithms can improve sea surface temperature (SST) accuracy.
Furthermore, [24] developed a practical split-window algorithm to estimate transmittance from water content for MODIS which requires transmittance and emissivity for LST retrieval. They concluded that their algorithm is not sensitive to water content and get higher accuracy if they can reasonably utilize the prior knowledge of water content. Standard atmospheric simulation and MODIS LST product are used to validate the algorithm and achieve the average accuracy of 0.32°C in LST retrieval for the case without error in both transmittance and emissivity estimations; and 0.37°C and 0.49°C, when the transmittance is computed from water by exponent fit and linear fit respectively. Compared with the MODIS LST product, the results from the analysis indicate that the algorithm is able to provide an accurate estimation of LST from MODIS data [24]. The corrections for atmospheric absorption, atmospheric emission, and surface emissivity are not easy to be defined, so the development of accurate LST algorithms is difficult [23,35]. The most robust algorithm for retrieving LST is the split window algorithm [24]. Unfortunately, the split-window algorithm cannot be applied for the analysis of Landsat 5 TM and Landsat 7 ETM+ data since it requires more than one thermal bands [22] (Otukei and Blaschke, 2012; [6,33]. In this regard, approaches based on single channels must be adopted [20,6]. In the Niger Delta, limited researches have been undertaken on LST measurement, for example [34] measured temperature at gas flaring sites in Rivers State using conventional method. [35] investigated geo-spatial dynamics of LST of PortHarcourt Metropolis and Environs using 1986, 2003 and 2018 Landsat data. They concluded that LST has increased from 30.2°C to 35.3°C with a variation of 5.1°C thereby leading to concentration of urban heat in some segments of the city. In addition, [36] worked on regional mapping of LST of South-South coastal settlements of Rivers State using Landsat data. They stated that LST decreased from the regional centre to the rural settlements ranging between 35.99 ºC to 22.12°C with a difference of 13.87°C. However, LST was observed scattering in different settlement spots especially in the Northern region showing that some settlements have started accelerating land surface modification activities. Furthermore, [37] examine the implications of land use change on LST and heat fluxes over three major cities (Benin City, Port Harcourt and Warri) in the Niger Delta region, using Landsat satellite data. The results showed a general increased in the mean LST, with an average of 1.43 °C increase between the year 2004 and 2015 over different urban land use. They stated that increase in LST appeared to be a result of changes in land use/land cover in the cities. The results further show that different land uses exhibit different degree of LST during both wet and dry seasons, with temperature nearly 2 °C higher in dry seasons compare to wet seasons. However, no paper has been published on the measurement of LST through satellite observations at gas flaring sites in the Niger Delta. Therefore, a valid way to bridge this gap and to have independent information on LST at flaring sites in the Niger Delta is from the use of Earth Observation (EO) satellites, owing to their continuous and repetitive observations, multispectral capabilities and synoptic coverage [38,39]. The research questions for this paper are: (1). How accurately can we retrieve LST from satellite based sensors in the Niger Delta? (2). What is the spatial variability and qualitative analysis of the LST retrieved from satellite based sensors in the Niger Delta? Based on these two research questions, this paper is aimed to discuss the ability of Landsat 5 TM and Landsat 7 ETM+ sensors to record LST at gas flaring sites in the Niger Delta. The specific objectives set out to answer the above research questions are: (1). Selection of appropriate flaring sites; (2). Estimation of emissivity values for each site land cover (LC); (3) Derivation of LST from atmospherically corrected Landsat 5 TM and Landsat 7 ETM+ data.
The implementation of methodologies based on the use of long-term time series of Landsat 5 TM and Landsat 7 ETM+ satellite data has demonstrated its potential and usefulness for this aim. In detail, satellite-based methods allow to map and monitor the distribution of LST over a large area, both in night-time [39,40,41] and daytime conditions [42,43]. The paper is organized as follows. Section 2 describes the study area and dataset used. In Section 3, the developed methodology was discussed; Landsat 5 TM and Landsat 7 ETM + data processing steps (Section 3.1), estimation of emissivity values for land covers (Section 3.2), retrieval of LST with atmospherically corrected Landsat 5 TM and Landsat 7 ETM + data (Section 3.3), characterization of spatial variability in LST (Section 3.4), and quantitative analysis of the LST (Section 3.5). In Section 4, the results on characterization of spatial variability in LST (Section 4.1), and that of the qualitative analysis of the LST (Section 4.2) are shown and discussed. Finally, Section 5 presents conclusions.

Study Area
In this study, four gas flaring sites in River States of the Niger Delta region, Nigeria were examined. The selection of the sites were based on the availability of Landsat data in the

Data Used
Six images of Landsat 5 TM and eleven images of Landsat 7 ETM+ dated January 17, 1986 to March 8, 2013 used were downloaded from the USGS Earth Resources Observation and Science (EROS) Data Centre website (http://earthexplorer.usgs.gov/) using the Glovis/Earth Explorer interface. The scenes with < 5 % contamination were selected for the study (Table 1); and all the sites are located within a single Landsat scene (Path 188, Row 057). Landsat 5 TM images are acquired in seven spectral bands while Landsat 7 ETM+is acquired in eight spectral bands and both have similar spatial resolution of 30 m for bands 1-5 and 7. Band 8 for Landsat 7 ETM+ is panchromatic with a spatial resolution of 15 m. The spatial resolution for Landsat 5 TM band 6 (thermal infrared) is 120 m while for Landsat 7 ETM+ is 60 m but both are resampled to 30 m pixels [20,45]. L1T is the processing level for all the scenes, which means systematic radiometric and geometric correction using ground control points (GCPs), and the digital elevation model (DEM) has been applied.
The problem of Scan Line Correction (SLC-off mode) with Landsat 7 sensor which started from 2004 that lead to loss of part of data in the scenes [46] was reduced to a minimum by setting one of the criteria for the selection of flare sites as the availability of data covering each facility. For these flaring sites, the size of the area investigated around the flare stacks with Landsat data is 12 km by 12 km, in order to include sufficient data for detail mapping of LST at each site so that processes not related to flaring could also be resolved.

Landsat 5 TM and Landsat 7 ETM+ Data Processing Steps
1) Verification of geo-location points using five ground control points that were selected over the Niger Delta using Google Earth. Two images each for both Landsat 5 TM and Landsat 7 ETM+ were uploaded into ArcGIS and the selected ground control points (GCPs) were identified. The comparison of coordinates of these controls obtained from both Google Earth and ArcGIS was carried out with a negligible difference found (1.0 × 10 -6 to 7.3 × 10 -6 m). This was taken as an acceptable error range for the geo-location of the imagery. 2) Removal of zero or out of range values from the data and their replacement with not a number (nan) to avoid divide by zero errors in calculations. MATLAB code was employed to process the data and to remove the zero and values at the upper and lower limits of the 8bit data range which cannot be distinguished from noise. 3) Radiometric calibration of the thermal band (6) Where: L λ = Spectral Radiance at the sensor's aperture in Wm⁻²sr⁻¹µm⁻¹; QCAL = the quantized calibrated pixel value in DN (Digital Number); LMIN λ = the spectral radiance that is scaled to QCALMIN in Wm⁻²sr⁻¹µm⁻¹; LMAX λ = the spectral radiance that is scaled to QCALMAX in Wm⁻²sr⁻¹µm⁻¹; QCALMIN = the minimum quantized calibrated pixel value (corresponding to LMIN λ ) in DN = 1 for LPGS (a processing software version) products; QCALMAX = the maximum quantized calibrated pixel value (corresponding to LMAX λ ) in DN = 255.
(4) Correction of the atmospheric effects for thermal band. The atmospheric correction parameters, upwelling radiance (L u ), downwelling radiance (L d ), and transmittance (τ) were obtained from Atmospheric Correction Parameters (ATMCORR) Calculator, a NASA web tool developed by [48]. Dark Object Subtraction (DOS) method [49,50] was adopted for the correction of atmospheric effects on the multispectral bands (1)(2)(3)(4) of Landsat data which is probably the most atmospheric correction method [51] reported in the literature.

Estimation of Emissivity ( ) Values
For this research, the four land cover (LC) types (vegetation, soil, built up and water) with their% identified at these sites during ground validation (July 02 to August 04, 2012 and February 04 to March 05, 2019) was derived when MATLAB code and cluster analysis method [52] were applied to the atmospherically corrected reflectance (bands 1-4) produced by Landsat 5 TM and Landsat 7 ETM+ [53]. The LC types was clarified using images held within Google Earth and Digital Global (http://browse.digitalglobe.com/imagefinder/public.do) such as WorldView-1 and 2, IKONOS pseudo-true colour images, Landsat imagery (bands 1-4 and 6), and Red, Green, Blue (RGB) pseudo-true colour composite images. The method used to estimate emissivity (Ɛ) value for LC types at these sites is based on the Ɛ value of four LC types present at each site. Each pixel LC types were considered for the entire site and their Ɛ values (both minimum (min) and maximum (max)) were taken from the literature. Mean of Ɛ values for LC types for each single pixel obtained from using their (min) and (max) values from the literature were calculated ( Table 2). Average of Ɛ values (min) and Ɛ (max) results were obtained for each pixel and the same procedure was repeated for all pixels in the selected 12 by 12 km area [54,20]. Therefore, the Ɛ value for each Landsat pixel is a combination of the Ɛ value of background features and that of any flare present within the pixel. The authors adopted an independent method of using LC types at each site for the correction of Ɛ value rather than Global Land Cover (GLC) data from USGS in order to ensure quality control primarily.

Retrieval of LST with Atmospherically Corrected Landsat 5 TM and Landsat 7 ETM+ Data
The theoretical basis for the LST measurement is Plank's radiation function, formulated as: Where: B (λ, T) is the spectral radiance of a blackbody in Wm -2 sr -1 µm -2 , and in practice, it is the emitted radiance of a ground object; λ is the wavelength in metres; T is temperature in Kelvin; C 1 = the first spectral constant =3.741775×10 -22 Wm 2 ; C 2 = the second spectral constant 1.4388 × 10 -2 mK; PI (π) = constant = 3.142 [6]. When B (λ, T) is measured, generally by a thermal sensor, the surface-leaving radiance (L λ ) or brightness temperature can be computed by inverting the Planck's radiance function as follows: L λ is the surface-leaving radiance or brightness temperature i.e. temperature corresponding to observed top of atmosphere (TOA) radiation for a black body, measured at TOA.
The approach for the calculation of LST, by first calculating L λ and substituting it into the Planck function and inverting the function to get the LST [20,65] was adopted for the study. The formula for computing L λ is: The L u , L d , and τwere applied to the calibrated at-sensor radiance band 6 (high gain) data to compute the L λ (Equation 4).
Then, the Planck equation was inverted using the calibration constants to derive LST (Equation 5).

Characterization of Spatial Variability in LST
This analysis is achieved through ArcGIS; and helps to fully characterise the 2D shape of each flare plume that enables a better understanding of the similarities and differences between individual plumes. It also helps to determine the best direction of the minimum and maximum LST slope for each site. For the figures, the pseudo-true colour images from the combination of Landsat bands 3, 2 and 1 as Red, Green and Blue (RGB) were used as the background map for the sites instead of the Google Earth images in order to avoid georeferencing errors that are associated with Google Earth images. The LST retrieved from Landsat 5 TM and Landsat ETM+ is uploaded into ArcGIS and processed to obtain the LST map in order to investigate its spatial variability for each site. In Figures 3-6, the bigger arrow and letter N at the upper corner of the right side shows the direction of the North, the small arrow pointed at the location of the flare within the site and Kelvin (K) in the legend is the International System Unit for temperature. LST was classified into 6 range group with a colour to represent each group. Pure red colour is for highest range of values, followed by light red, light brown, deep orange, light orange and yellow colours respectively.

Qualitative Analysis of the LST
The primary aim of this analysis is to distinguish the flare stack position from other parts of the facility within the site using the retrieved LST. Some of the factors that can affect LST at these flaring sites include the rate of gas burning, stack height, facility size, vegetation type, vegetation density and time of observation (month) i.e. rainy season or dry season. Those that are available for this study are facility size, stack height and time. Each flaring site is made up of 400 by 400 Landsat 5 and Landsat 7 resampled pixels, with the flare stack positioned at the centre of the site. The transect plots of LST was carried out by plotting at the pixel 200 from the South through the centre of the flare stack to the North at 200 pixels (Figure 7). South-North is the direction of prevailing wind in the Niger Delta [20]. For all these facilities, the pixels for the flare stacks were found to have higher LST values than that of the surrounding pixels. In addition, in some cases LST was elevated in adjoining pixels, suggesting either a warming effect of the flare on surrounding structures or reflecting structures. In a few instances, such as Bonny LNG, some thermal band pixels are also brighter possibly as a result of other facilities such as metal oil storage tanks that absorb heat from both the flare and sun and zero flare at the time of satellite overpass. The author clarified this situation using a combination of images from Google Earth and results from Landsat visible bands (1)(2)(3)(4) to ensure that such thermal band pixels were not wrongly interpreted as the flaring source.

Results and Discussion
In this section, LST maps and LST transect plots (South-North direction) obtained from Landsat 5 TM and Landsat 7 ETM+ for each site are presented in Figures 3-7. Figures 3A-6A, show the 2D plot for LST for four flaring sites examined. Figures 3B-6B show the 2D plot for LST in 6 different layers while Figures 3C-6C are Figures 3B-6B with additional contours (contour interval of 0.5 K) in order to show the extent and nature of variation in LST within each site. For all four sites investigated, the spatial analysis of LST shows that the flare sources gives the highest LST, followed by the next adjoining pixels surrounding the flare and continue in that order.  For Eleme Refinery II (Figures 3 A, B and C), six classes of LST values obtained are 316-320 K for pure red points, 312-316 K for light red points, 308-312 K for light brown points, 304-308 K for deep orange points, 300-304 K for light orange points and 296-300 K for yellow points ( Table  3). The plot shows that the flare source has the highest range of values of LST (pure red points) that spread from the centre of the plume downward and toward South-West, West, North-East and East directions. LST recorded in the South-East direction is a mixture of the second (light red points), third (light brown points) and sixth (yellow points) classes of LST with a few of red points (Figure 3 A)  Six classes of the range of LST values for Bonny LNG (Figures 4 A, B and C) are 332-358 K (pure red points), 320-332 K (light red points), 317-320 K (light brown points) and 313-317 K (deep orange points), 309-313 K (light orange points) and 301-309 K (yellow points) ( Table 3) (Table 3).   Figures 6 A, B and C as 312-319 K (pure red points), 306-312 K (light red points), 303-306 K (light brown points), 300-303 K (deep orange points), 296-300 K (light orange points) and 290-296 K (yellow points) ( Table  3). The result in Figures 6 A, B and C shows the extent of the plume from the flare stack and its immediate surrounding pixels gives the highest range of LST values for the scene. The entire result is dominated by LST values from the fourth and fifth classes. Also, there are up to 6 LST values from the third class in the North direction of the scene (Figure 6 A). The result suggests that the South direction of the prevailing wind for the Niger Delta could not have significant effect on the flare. The dimension of the plume within the site for this scene is 9 by 11 pixels.

Characterization of Spatial Variability in LST
In summary, average of LST values recorded for 6 different layers for each flaring site is presented in Table 4.
Bonny LNG recorded the highest LST value (345.0 K), followed by Obigbo Flow Station (333.5 K), Eleme Petroleum Company Refinery II (318.0 K) and Umudioga Flow Station (315.5 K) respectively. It is observed that the size and shape of the plume differ from one oil facility to another. For example, the sizes of the plume for the bigger facilities are 11 by 13 pixels for Eleme Refinery II, 16 by 15 pixels for Bonny LNG respectively. For small facilities the sizes of the plume recorded are 9 by 7 pixels for Obigbo Flow Station, 9 by 11 pixels for Umudioga Flow Station. These results show that the sizes of the plumes from the bigger oil facilities are larger than those obtained from small facilities. Therefore, the result suggests that the major factors that determine the size of the plume are the size of facility, volume of burning gas and the rate of its burning at the time of satellite overpass.  Tables 3 and 4); and followed by the immediate surrounding pixels in the order of their closeness to the flare source. The range of values of LST recorded for each site is different; and it is observed from the results that the size and shape of the plume also differ from one site to another. The sizes of the plumes from the bigger oil facilities are larger than those obtained from and small facilities. However, the range of values of LST and the average of LST recorded does not depend on the largeness of the oil facility within the site. For example, a small Obigbo Flow Station (650 × 650) m produced average value of LST of 333.5 K which is higher than that of bigger facility such as Eleme Refinery II (2.2 × 1.3) with 318 K as the average value of LST. Figure 7 show the results of the South-North transect plots of the LST recorded for each site. Similarly, all the flare stacks produced LST of highest values which confirm those points as the hottest spot within the site as obtained by the spatial analysis discussed in section 4.1.

Qualitative Analysis of LST
In summary, Landsat 5 TM and Landsat 7 ETM+ can retrieve LST accurately in the Niger Delta. The results presented in Figures 3-6 for characterization of variability in LST and Figure 7 for qualitative analysis of LST shows that each flare source produced the highest LST, followed by the adjoining pixels and continue in that order. The results presented in Tables 3 and 4 shows that Obigbo Flow Station which is a smaller facility recorded a higher range of LST (330-337) K and average of LST (333.5) K compared to Eleme II (316-320) K and (318.0) K which is a larger facilities. This suggests that the volume of gas flared and the rate at which the gas is flared at the time of satellite overpass determined the amount of LST recorded by the satellite not the size of the facility. Also, the larger the size of the facility, the larger the size of the plume recorded by the satellite. Furthermore, the spatial variability of LST for all sites varies, it differ from one site to another. This suggests that apart from difference in the size of facilities, the volume of gas flared and the rate of gas flaring varies for all sites. Since it can be confirmed from the results of this research that Landsat 5 TM and Landsat 7 ETM+ can accurately retrieved LST at flaring sites in the Niger Delta, the following recommendations are made: The primary restriction to this study is the lack of satellite derived and in-situ data. Data from the two satellites owned by the Nigerian Government could not be accessed. Therefore, Nigerian Government should ensure that these data are available to educational and research institutions that will need them for research.
Lack of direct and open access for the measurement of insitu data at oil and gas facility sites in the Niger Delta because of the security issues; lack of data on the volume of gas flare and the rate of flaring gas from these sites are further limiting factors in this study. Therefore, the government should make a provision for the policies that enforce multi-national oil companies to declare information on their oil and gas exploration and exploitation activities to the general public especially to stakeholders and organizations involved in the oil and gas business.
The provision of an enabling environment and sufficient funding for scientific research on oil and gas related disciplines such as gas flaring and to fully assess it impacts on vegetation, biodiversity and ecosystem, and ensure ways of mitigation. Staffs training especially for those that are in charge of implementation of policy relating to oil and gas production processes and land economics.

Conclusion
Based on the results from Figures 3-7 and Tables 3 and 4, Landsat 5 TM and Landsat 7 ETM+ sensors could be used to record LST for the flaring sites in the Niger Delta. LST retrieved from both sensors for the flare hotspots are the highest values compared to other locations within the processing sites, which was clearly shown through GIS spatial analysis and the transects plots. Also, the closer is the distance to the flare, the higher is the temperature and vice versa. Based on these results, it can be concluded that satellite based sensors, such as Landsat TM and ETM+, have the ability to record LST at gas flaring sites in the Niger Delta, Nigeria.