Atmospheric Correction of Thermal Infrared Landsat Images Using High-Resolution Vertical Proﬁles Simulated by WRF Model †

: Atmospheric proﬁles are key inputs in correcting the atmospheric effects of thermal infrared (TIR) remote sensing data for estimating land surface temperature (LST). This study is a ﬁrst insight into the feasibility of using the weather research and forecasting (WRF) model to provide high-resolution vertical proﬁles for LST retrieval. WRF numerical simulations were performed to downscaling NCEP climate forecast system version 2 (CFSv2) reanalysis proﬁles, using two nested grids with horizontal resolutions of 12 km (G12) and 3 km (G03). We investigated the use of these proﬁles in the atmospheric correction of TIR data applying the radiative transfer equation (RTE) inversion single-channel approach. The MODerate resolution atmospheric TRANsmission (MODTRAN) model and Landsat 8 TIRS10 (10.6–11.2 µ m) band were taken for the method application. The accuracy evaluation was performed using in situ radiosondes in Southern Brazil. We included in the comparative analysis the atmospheric correction parameter calculator (ACPC; NASA) web tool and proﬁles directly from the NCEP CFSv2 reanalysis. The atmospheric correction parameters from ACPC, followed by CFSv2, had better agreement with the ones calculated using in situ radiosondes. When applied into the RTE to retrieve LST, the best results (RMSE) were, in descending order: K), and WRF G03 (0.82 K). The ﬁndings suggest that increasing the horizontal resolution of reanalysis proﬁles does not particularly improve the accuracy of RTE-based LST retrieval. However, the WRF results are yet satisfactory and promising, encouraging further assessments. We endorse the use of the well-known ACPC and also recommend the NCEP CFSv2 reanalysis proﬁles for TIR remote sensing atmospheric correction and LST single-channel retrieval.


Introduction
Land surface temperature (LST) is a key parameter in a wide variety of environmental applications. It is closely connected and plays an important role in Earth's surfaceatmosphere interactions at both local and global scales [1][2][3][4]. Thermal infrared (TIR) remote sensing is an outstanding way of obtaining the LST at regional and global scales [5][6][7]. Nevertheless, the spectral radiation measured by TIR sensors onboard satellites is determined not only by the surface parameters (emissivity and temperature) but also by the atmosphere effects (mainly due to water vapor) [6,8]. These atmospheric perturbations The Porto Alegre International Airport (SBPA), Rio Grande do Sul State, Brazil was selected as the study area. The SBPA includes a radiosonde station, which makes this site useful for studies that aim to validate atmospheric profiles. The selected area covers the official limits of the Anchieta district, with an area of around 9.2 km 2 ( Figure 1). The SBPA station is located at 30.00° S and 51.18° W, 3.0 m above mean sea level. Radiosondes are launched twice a day, at 00:00 and 12:00 UTC. We used the 12:00 UTC radiosonde profiles, as this is the closest to the Landsat 8 crossing time over the study site (~13 UTC). This dataset allows characterizing the vertical structure of the atmosphere with profiles of air temperature, pressure, humidity with up to 99 vertical levels. The radiosonde observations, as well as the parameters calculated from them, were used as the ground truth for the assessments.

Satellite and Reanalysis Data
Landsat 8 carries a two-sensor payload: the operational land imager (OLI) that has nine reflective bands and the thermal infrared sensor (TIRS) with two bands in the TIR region. We acquired all Landsat 8 images (Collection 1) available under daily clear-sky conditions over the study site (Path-Row 221-81) from 2013 to 2019, which resulted in a total of 27 scenes, henceforward refer as case days 1-27. The full swath Landsat data were subset to a 10,184 pixels region covering the study area of Figure 1. The TIRS band 10 (10.60-11.19 µm) is used for RTE-based LST retrieval and OLI bands 4-red (0.64-0.67 µm) and 5-near-infrared (NIR) (0.85-0.88 µm) for land surface emissivity estimation.
Additionally, we used the NCEP climate forecast system version 2 (CFSv2) [32] reanalysis data from the 6-hourly product as initial and boundary conditions for the WRF simulations. NCEP CFSv2 reanalysis data is produced using the NCEP global forecasting system (GFS) atmospheric model and the grid point statistical interpolation (GSI) analysis system with three-dimensional variational data assimilation (3D-Var). It is arranged in grids with a horizontal resolution of 0.5° × 0.5° in 37 vertical (pressure) levels (1000-1 mbar), and of 0.205° for surface parameters. Besides, profiles retrieved directly from NCEP CFSv2 were included in the analysis, to assess the WRF model horizontal resolution downscaling performance. These profiles were extracted from the grid point closest to the SBPA station.

Satellite and Reanalysis Data
Landsat 8 carries a two-sensor payload: the operational land imager (OLI) that has nine reflective bands and the thermal infrared sensor (TIRS) with two bands in the TIR region. We acquired all Landsat 8 images (Collection 1) available under daily clear-sky conditions over the study site (Path-Row 221-81) from 2013 to 2019, which resulted in a total of 27 scenes, henceforward refer as case days 1-27. The full swath Landsat data were subset to a 10,184 pixels region covering the study area of Figure 1. The TIRS band 10 (10.60-11.19 µm) is used for RTE-based LST retrieval and OLI bands 4-red (0.64-0.67 µm) and 5-near-infrared (NIR) (0.85-0.88 µm) for land surface emissivity estimation.
Additionally, we used the NCEP climate forecast system version 2 (CFSv2) [32] reanalysis data from the 6-hourly product as initial and boundary conditions for the WRF simulations. NCEP CFSv2 reanalysis data is produced using the NCEP global forecasting system (GFS) atmospheric model and the grid point statistical interpolation (GSI) analysis system with three-dimensional variational data assimilation (3D-Var). It is arranged in grids with a horizontal resolution of 0.5 • × 0.5 • in 37 vertical (pressure) levels (1000-1 mbar), and of 0.205 • for surface parameters. Besides, profiles retrieved directly from NCEP CFSv2 were included in the analysis, to assess the WRF model horizontal resolution downscaling performance. These profiles were extracted from the grid point closest to the SBPA station.

WRF Model Configuration
The WRF Model version 4.1.2 with the advanced research WRF (ARW) dynamical solver [27] was used to perform high-resolution numerical simulations. We configured the WRF domains with two nested grids, in one-way mode, centered at SBPA station, with horizontal resolutions of 12 km (G12) and 3 km (G03) (4:1 parent grid ratio), and 33 sigma vertical levels with 50 hPa top pressure value. The WRF configurations used are summarized in Table 1 based on Santos and Nascimento (2016) [33] for the same station. The model was run in the below settings for each case day. The profiles resulted from simulations were retrieved from the model grid point closest to the SBPA station We conduct simulations of 24 h duration and the results were extracted at 12:00 UTC, to match with the local radiosonde observations. Thus, the first 12 h of the simulation was considered for a spin up time.

Land Surface Emissivity Estimation
Land surface emissivity (LSE) is one of the key parameters to retrieve LST from remote sensing data [34,35]. Among the methods for LSE retrieval from space, normalized difference vegetation index (NDVI)-based ones are operational and the most applied, with satisfactory results [36][37][38][39]. Bonafoni (2020a, 2020b) [40,41] examined the influence of six NDVI-based LSE models on the performance of LST retrieval. Based on their results, we used the NDVI threshold method (NDVI THM ) proposed by Sobrino et al. (2008) [36].
The NDVI is obtained from Equation (1): where ρ NIR and and ρ R are the reflectances of NIR and red bands, respectively. From NDVI, it is possible to calculate the fractional vegetation cover (P V ) from Equation (2) [42]. The P V is an important factor in the LSE estimation.
The NDVI THM estimates LSE considering three different cases as presented in Equation (3), for Landsat 8 [40]. For NDVI < 0.2, the pixel is considered as bare soil, and the emissivity is calculated using the reflectance of the red band. In the second case (0.2 ≤ NDVI ≤ 0.5), the pixel is considered composed of a mixture of bare soil and vegetation and the LSE depends on the P V value. The pixels with NDVI values higher than 0.5 are considered as fully vegetated areas and the emissivity is assumed as 0.99.
where ε is the Land Surface Emissivity (LSE).

Atmospheric Parameters Calculation with MODTRAN and ACPC
The present study used the MODTRAN (MODerate resolution atmospheric TRANsmission) 4 v3r1 [43] to estimate the three atmospheric correction parameters (i.e., atmospheric transmittance, upwelling atmospheric radiance, and downwelling atmospheric radiance) in the Landsat TIR spectrum. We introduced into the MODTRAN as input vertical profiles of pressure, air temperature, and relative humidity from: (i) SBPA radiosonde; (ii) NCEP CFSv2 reanalysis; (iii) WRF G12; and (iv) WRF G03.
The methodology of Barsi et al. (2003) [12], to fulfill the profiles, was adopted. To predict space-reaching atmospheric parameters, the MODTRAN requires atmospheric profiles reaching "space", or 100 km above sea level. Since the radiosondes and NCEP CFSv2 are from the surface to about 30 and 50 km, respectively, the upper atmosphere layers (to 100 km) are extracted from the MODTRAN standard atmospheres and pasted onto our site-specific profiles. We take on the standard mid-latitude summer profile for case days in hot seasons (spring and summer) and the mid-latitude winter profile for those in cold seasons (autumn and winter). This results in surface-to-space vertical profiles of air temperature, pressure, and water vapor. These completed profiles are those that are inserted into a MODTRAN input file and then processed [10,12].
Moreover, we include in the comparative analysis the atmospheric parameters estimated by the well-known atmospheric correction parameter calculator (ACPC; NASA) web tool [10,12]. The ACPC uses NCEP reanalysis profiles (1 • × 1 • horizontal resolution and 28 vertical levels), MODTRAN code, and a suite of integration algorithms to directly provide the AC parameters for a particular date, time, and location inputted. The option of using the atmospheric profile from the closest integer coordinate to the inputted location (SBPA station) was set. The mid-latitude standard upper profiles varied according to the season of each case day.

Radiative Transfer Equation (RTE) Based LST Retrieval Method
The inverse solution of RTE [11] is a direct and a priori the most appropriate procedure for LST retrieval using a single TIR band [44]. The RTE applied to a particular TIR band/wavelength (λ) can be simplified and given by: where L sen λ (W·m −2 ·sr −1 ·µm −1 ) is the at-sensor (TOA) spectral radiance of the corresponding TIR band (in this paper, TIRS10), ε λ refers to the LSE (dimensionless), B λ (W·m −2 ·sr −1 ·µm −1 ) is the blackbody radiance, T s (Kelvin) represents the LST, L ↓ λ , and L ↑ λ (W·m −2 ·sr −1 ·µm −1 ) refer to the downwelling and upwelling radiances, respectively, and τ λ is the atmospheric transmittance (dimensionless). So, the emitted radiance for a black body at a temperature T s is given by the inversion of Equation (4): the T s is calculated by inverting Planck's law in Equation (5) and the LST (T s ) from Landsat 8 TIRS10 is estimated as: where K 1 and K 2 refers to calibration constants, whose values for the Landsat 8 TIRS10 are 774.89 W·m −2 ·sr −1 ·µm −1 and 1321.08 K, respectively [45]. Henceforth, the spectral notation (λ) will be omitted, since here only a single TIR band is used.

Metrics for Performance Evaluation
To evaluate the performance of the WRF model and other profiles we take into account the atmospheric parameters and LST images. The SBPA radiosondes are the available in situ observations. Hence, the AC parameters and LSTs calculated using SBPA profiles are considered our references. To perform the comparative assessment, Pearson's correlation coefficient (R), bias (mean error), mean absolute error (MAE), and root mean square error (RMSE) were used as statistical criteria. These metrics are widely employed to evaluate and compare models.

Evaluation of Atmospheric Parameters
The AC parameters (τ, L ↑ , and L ↓ ) calculated with the different sources of estimated vertical profiles (CFSv2, WRF G12, WRF G03, and ACPC) were compared against those using observational SBPA radiosondes. In Table 2, the accuracy of atmospheric parameters estimations is presented. All the profile sources provide AC parameters with high correlation coefficients (all greater than 0.9) concerning the reference (SBPA). The R values of ACPC, followed by CFSv2, are slightly better. There is a general but small tendency of overestimating the transmittance values. On the other hand, the atmospheric radiances tend to be underestimated, except for ACPC downwelling. The smallest biases were from the WRF for all three parameters. The largest ones were from CFSv2. Concerning MAE and RMSE, the best results were from ACPC followed by CFSv2. The higher negative bias of the atmospheric radiances with NCEP CFSv2 in our finds may be due to that these reanalysis profiles have the lowest level at 1000 hPa, which means around 60-250 m for SBPA station in analyzed case days. So the lowest layer of the atmosphere (which typically presents the largest water vapor content and warmest temperature) is neglected in these profiles [5,14]. We had tried to reduce this limitation by downscaling the CFSv2 reanalysis profiles with the WRF model. The WRF profiles bring the first level at around 1 m above the surface.
Our RMSE range is in agreement with those found for reanalysis profiles (eight different products) in Meng and Cheng (2018) [5] for the three parameters. The authors analyzed 30,715 atmospheric profiles in 163 stations around the globe. The RMSE values of Yang et al. (2020) [46] were overall slightly lower than ours. Their assessment included seven profiles sources over Europe.
The errors in the estimation of atmospheric parameters for each case day are shown in Figure 2. It shows that despite ACPC presented the best overall metrics, none of the profile sources outperforms in all cases. For instance, in 8 of the 27 case days, one of the WRF profiles had the best results in calculating the downwelling radiance. In case day 23, the WRF model produced the largest errors, whereas in cases such as 26 the model successfully reduced the largest error of the driving reanalysis data.
The higher negative bias of the atmospheric radiances with NCEP CFSv2 in our finds may be due to that these reanalysis profiles have the lowest level at 1000 hPa, which means around 60-250 m for SBPA station in analyzed case days. So the lowest layer of the atmosphere (which typically presents the largest water vapor content and warmest temperature) is neglected in these profiles [5,14]. We had tried to reduce this limitation by downscaling the CFSv2 reanalysis profiles with the WRF model. The WRF profiles bring the first level at around 1 m above the surface.
Our RMSE range is in agreement with those found for reanalysis profiles (eight different products) in Meng and Cheng (2018) [5] for the three parameters. The authors analyzed 30,715 atmospheric profiles in 163 stations around the globe. The RMSE values of Yang et al. (2020) [46] were overall slightly lower than ours. Their assessment included seven profiles sources over Europe.
The errors in the estimation of atmospheric parameters for each case day are shown in Figure 2. It shows that despite ACPC presented the best overall metrics, none of the profile sources outperforms in all cases. For instance, in 8 of the 27 case days, one of the WRF profiles had the best results in calculating the downwelling radiance. In case day 23, the WRF model produced the largest errors, whereas in cases such as 26 the model successfully reduced the largest error of the driving reanalysis data.  On the whole, ACPC obtained better results than the other profiles for all three atmospheric parameters. Then, comes the CFSv2 reanalysis. Table 2 points out that no significant statistical differences were found between the parameters' accuracies from the WRF grids G12 and G03. It suggests that computation costs can be saved by using profiles from a WRF domain with coarse horizontal resolution. Despite other aims, scholars have already reported this kind of result with WRF model grids [24,[47][48][49][50].

Application to RTE-Based LST Retrieval
To further assess the different atmospheric profiles, the retrieved LSTs by RTE inversion with atmospheric parameters, Landsat 8 TIRS10 radiance, and NDVI THM emissivity were intercompared. Once more, the LST images that used SBPA profiles were assumed as reference data. Except for the atmospheric parameters calculated from the different profile sources, the other variables in Equation (6) were the same for each pixel of the scenes. Therefore, the differences in LST values are due to the discrepancies among the profiles [46].
Histograms of LST errors take in all the 10,184 pixels in the study area of the 27 case days are shown in Figure 3. These represent the frequency distribution of the errors in the retrieval of LST using the different atmospheric profile sources. For all profiles, more than 50% of the errors are of ±1 K. Yang et al. (2020) [46] also found most of the LST differences in this range, using seven different reanalysis and satellite-derived profiles. The histograms in Figures 3b and 4c indicate that WRF profiles tend to overestimate the LST. Whereas ACPC tends to underestimate it (Figure 3d). Using WRF profiles, LST errors can reach more than 4 K, although in a very small number of cases. For ACPC and both WRF grids, the error range that occurs most often is between 0 and −1 K. The distribution of CFSv2 LST errors is more symmetrical than in the other profiles. Table 3 summarizes the metrics of the LST retrieval comparative analysis. Overall, the LST values of the four profile sources analyzed in this study were found in good agreement with the reference. All of them showed a very strong correlation and relatively low bias, MAE, and RMSE values. CFSv2, WRF G12, and WRF G03 presented an average positive bias and ACPC a negative one. This corroborates with the histogram analysis in Figure 3. The mean error criteria (MAE and RMSE) indicate that the profiles with the best performance in the RTE-based LST retrieval are, in descending order: CSFv2, ACPC, WRF G12, and WRF G03. The differences between CFSv2 and ACPC overall MAE and RMSE values were very small. The same for WRF G12 and G03.  On the whole, ACPC obtained better results than the other profiles for all three atmospheric parameters. Then, comes the CFSv2 reanalysis. Table 2 points out that no significant statistical differences were found between the parameters' accuracies from the WRF grids G12 and G03. It suggests that computation costs can be saved by using profiles from a WRF domain with coarse horizontal resolution. Despite other aims, scholars have already reported this kind of result with WRF model grids [24,[47][48][49][50].

Application to RTE-Based LST Retrieval
To further assess the different atmospheric profiles, the retrieved LSTs by RTE inversion with atmospheric parameters, Landsat 8 TIRS10 radiance, and NDVI THM emissivity were intercompared. Once more, the LST images that used SBPA profiles were assumed as reference data. Except for the atmospheric parameters calculated from the different profile sources, the other variables in Equation (6) were the same for each pixel of the scenes. Therefore, the differences in LST values are due to the discrepancies among the profiles [46].
Histograms of LST errors take in all the 10,184 pixels in the study area of the 27 case days are shown in Figure 3. These represent the frequency distribution of the errors in the retrieval of LST using the different atmospheric profile sources. For all profiles, more than 50% of the errors are of ±1 K. Yang et al. (2020) [46] also found most of the LST differences in this range, using seven different reanalysis and satellite-derived profiles. The histograms in Figures 3b and 4c indicate that WRF profiles tend to overestimate the LST. Whereas ACPC tends to underestimate it (Figure 3d). Using WRF profiles, LST errors can reach more than 4 K, although in a very small number of cases. For ACPC and both WRF grids, the error range that occurs most often is between 0 and −1 K. The distribution of CFSv2 LST errors is more symmetrical than in the other profiles.    Table 3 summarizes the metrics of the LST retrieval comparative analysis. Overall, the LST values of the four profile sources analyzed in this study were found in good agreement with the reference. All of them showed a very strong correlation and relatively low bias, MAE, and RMSE values. CFSv2, WRF G12, and WRF G03 presented an average positive bias and ACPC a negative one. This corroborates with the histogram analysis in Figure 3. The mean error criteria (MAE and RMSE) indicate that the profiles with the best performance in the RTE-based LST retrieval are, in descending order: CSFv2, ACPC, WRF G12, and WRF G03. The differences between CFSv2 and ACPC overall MAE and RMSE values were very small. The same for WRF G12 and G03. Comparing with previous studies that evaluated the application of different atmospheric profiles in the LST retrieval, Meng and Cheng (2018) [5] reported overall LST RMSE values higher than ours using the Landsat 8 TIRS10 band and eight different reanalysis profile sources analyzed around the globe. All their eight average RMSEs were larger than 1 K. They also indicated an average tendency to overestimate the LST. In a study by Yang et al. (2020) [46], RMSEs ranged between 0.6-1 K comparing LST retrieved from Landsat 8 TIRS10 over Europe. Calculating LST from three MODIS thermal bands, Pérez-Planells et al. (2015) [51] showed RMSEs between 0.6 and 0.9 K using ACPC/NCEP and between  Comparing with previous studies that evaluated the application of different atmospheric profiles in the LST retrieval, Meng and Cheng (2018) [5] reported overall LST RMSE values higher than ours using the Landsat 8 TIRS10 band and eight different reanalysis profile sources analyzed around the globe. All their eight average RMSEs were larger than 1 K. They also indicated an average tendency to overestimate the LST. In a study by Yang et al. (2020) [46], RMSEs ranged between 0.6-1 K comparing LST retrieved from Landsat 8 TIRS10 over Europe. Calculating LST from three MODIS thermal bands, Pérez-Planells et al. (2015) [51] showed RMSEs between 0.6 and 0.9 K using ACPC/NCEP and between 1.3 and 3 K for MOD07 profiles, depending on the band and the altitude of the study sites in Spain. Figure 4 displays the distribution of LST bias and RMSE along the case days. Figure 4a evidences the ACPC and WRF settings average tendency of under and overestimate the LST, respectively. We found better RMSE values using CFSv2 in 11 of the 27 case days, and in 8 using ACPC. Although the metrics of mean errors are greater for WRF profiles, in 11 case days the model improves the results of the reanalysis. Besides, the finer grid (G03) succeeded in downscaling the G12 11 times. The fact is that the largest errors were achieved when using the WRF model (e.g., case days 23, 19, and 17), and it contributes to the higher overall RMSE values. In general, the days with larger errors in atmospheric parameters ( Figure 2) are the days with larger LST RMSEs, as evinces case 23. Conversely, in case 22, the errors in atmospheric parameters using ACPC are less than those using WRF profiles, but the highest LST RMSE on this day is with ACPC. Jiménez-Muñoz et al. (2010) [8] advocate that cases like this may be explained due to compensation among AC parameters errors; for instance, a significant positive difference in transmittance and significant but negative differences in the atmospheric radiances.
In summary, the essay of downscaling the horizontal resolution of reanalysis data from 0.5 • (~56 km) to 12 km and so to 3 km, aiming to reduce errors in the calculation of atmospheric parameters and hence in LST retrieval, did not perform as theoretically expected. With WRF simulations, we also improved the vertical resolution in the lowest atmospheric levels. However, no significant improvement was found using the WRF profiles in the AC. In some cases, using a finer grid resolution profile resulted in greater uncertainties in atmospheric parameters and LST estimation. Rosas et al. (2017) [7] reported that the higher vertical resolution of NCEP and ECMWF profiles in their study did not seem to play a significant role in the atmospheric correction. Even if, naturally, data with higher resolution tend to better represent the atmosphere parameters, it is not a strictly direct relationship [5,46]. Furthermore, the ACPC that uses profiles with 1 • × 1 • (~100 km) horizontal resolution showed good results. Although these profiles have a coarser horizontal resolution, previous studies have been finding satisfactory results using the ACPC, even surpassing other methods [8,16,38,52,53]. It is important to note that in this study the WRF resulting profiles were extracted at 12 UTC to match with the available radiosonde data. Nevertheless, it could be set for the exact time of the satellite overpass. In ACPC and, in general, for reanalysis profiles, this time synchronization is done through linear interpolation. Which may not be the most appropriate strategy for sampling weather fronts and diurnal heating cycles [5,10,12].

Conclusions
Vertical atmospheric profiles are key inputs in the atmospheric correction for estimating LST using the RTE inversion single-channel approach. This study evaluated the use of the WRF numerical model to simulate high-resolution profiles improving horizontal, temporal, and vertical resolutions of NCEP CFSv2 reanalysis data. The obtained results showed that the ACPC provided atmospheric parameters in better agreement with those calculated using radiosondes. The second-lowest differences were using CSFv2 profiles. No significant statistical differences were found between the parameters from the two WRF grids. None of the profile sources outperformed in all case days analyzed. The overall metrics of WRF profiles were influenced by some cases with large errors. Concerning retrieved LST values, CFSv2 had the best results. With an RMSE of 0.55 K, it was slightly more accurate than ACPC (RMSE of 0.56 K). WRF G12 and G03 showed RMSE values of 0.79 and 0.82 K, respectively. On balance, all the profile sources presented relatively good results in estimating the LST.
From the above findings, our main conclusion is that there is no special need to increase the horizontal resolution of reanalysis profiles aiming at general RTE-based LST retrieval. We recommend the use of NCEP CFSv2 profiles for these applications. Moreover, our results reinforce the ACPC validity and feasibility, which is free of charge. Even though the overall statistical metrics for WRF profiles were inferior, their results were satisfactory, both in the estimation of atmospheric parameters and LST values. Despite some studies used the WRF model to simulate the skin temperature [54,55], to the best of our knowledge, this paper is the first effort applying the WRF model to aid the atmospheric correction of thermal remote sensing data. Its use showed potential, and our findings encourage further validation. Our proposal joins the background for studies combining TIR satellite images and high-resolution NWP models.