Land Surface Temperature Retrieval Using High-Resolution Vertical Profiles Simulated by WRF Model ^†

Abstract

This work gives a first insight into the potential of the Weather Research and Forecasting (WRF) model to provide high-resolution vertical profiles for land surface temperature (LST) retrieval from thermal infrared (TIR) remote sensing. WRF numerical simulations were conducted to downscale NCEP Climate Forecast System Version 2 (CFSv2) reanalysis profiles, using two nested grids with horizontal resolutions of 12 km (G12) and 3 km (G03). We investigated the utility of these profiles for the atmospheric correction of TIR data and LST estimation, using the moderate resolution atmospheric transmission (MODTRAN) model and the Landsat 8 TIRS10 band. The accuracy evaluation was performed using 27 clear-sky cases over a radiosonde station in Southern Brazil. We included in the comparative analysis NASA’s Atmospheric Correction Parameter Calculator (ACPC) web-tool and profiles obtained directly from the NCEP CFSv2 reanalysis. The atmospheric parameters from ACPC, followed by those from CFSv2, were in better agreement with parameters calculated using in situ radiosondes. When applied into the radiative transfer equation (RTE) to retrieve LST, the best results (RMSE) were, in descending order: CFSv2 (0.55 K), ACPC (0.56 K), WRF G12 (0.79 K), and WRF G03 (0.82 K). Our findings suggest that there is no special need to increase the horizontal resolution of reanalysis profiles aiming at RTE-based LST retrieval. However, the WRF results were still satisfactory and promising, encouraging further assessments. We endorse the use of the well-known ACPC and recommend the NCEP CFSv2 profiles for TIR atmospheric correction and LST single-channel retrieval.

1. Introduction

Land surface temperature (LST) is one of the essential climate variables (ECVs) of the Global Climate Observing System (GCOS) [1,2]. It is closely connected to Earth–atmosphere interactions, playing a pivotal role in surface energy and water balances at both local and global scales [3,4,5]. Therefore, LST is a key parameter in a wide range of environmental applications [6]: urban heat island studies [7]; numerical weather prediction [8]; agricultural, forests and drought monitoring [9,10,11]; monitoring of geothermal activity and natural hazards [12,13]; evapotranspiration estimation [14,15]; fire detection [16,17]; water resource management [18,19]. It is worth mentioning that the LST should not be confused with the near-surface air temperature (typically measured by meteorological/in situ stations); the LST refers to the so-called “skin temperature” [10].

Thermal infrared (TIR) remote sensing is a one-off way of obtaining the LST at regional and global scales [20,21,22]. However, the spectral radiance measured by the TIR sensors on board satellites is influenced not only by the surface parameters (emissivity and temperature) but also by the composition and structure of the atmosphere (mainly water vapor) [23,24]. Thus, these atmospheric effects must be removed for the appropriate use of TIR remote sensing data in temperature research applications [4,25,26]. The atmospheric correction (AC) is, in general terms, the conversion of top-of-the-atmosphere (TOA) to ground-level measurements. Neglecting the AC leads to systematic errors in the LST estimation for any atmosphere [24,27].

One of the widely applied methodologies for AC and LST retrieval is the physics-based radiative transfer equation (RTE) method [28]. It involves a simple inversion of the RTE for a particular channel and can provide theoretically accurate LST retrieval [23]. The RTE approach requires vertical atmospheric profiles (air temperature, water vapor, and pressure). This information is introduced into a radiative transfer model (RTM) to calculate the three atmospheric parameters necessary for AC: atmospheric transmittance, upwelling atmospheric radiance, and downwelling atmospheric radiance [24,29]. In situ radiosonde profiles launched simultaneously with the satellite overpass are ideal for AC [30,31]. Nevertheless, this kind of profile is unavailable under most realistic conditions [20,32,33], as local radiosonde launching has a significant financial cost [23,24]. Radiosonde stations with atmospheric-profile databases are launched daily around the globe (e.g., in airports). However, these data may be either non-synchronous or too far away from the scene footprint [4,34]. Hence, local radiosondes are mainly suitable for particular local studies and validation at specific sites [21,35].

Therefore, atmospheric profile products from different sources have been used as surrogate radiosondes in TIR atmospheric correction, resulting in LSTs with acceptable accuracy [4,26,36,37,38,39]. Satellite-derived profiles overcome the spatial limitations of local radiosondes. These are available at the satellite pixel scale, providing data over a large spatial extent. Although they have a high spatial resolution, the satellite-derived profiles can be compromised by low temporal coverage (only at the satellite overpassing time) [21,40]. Thus, the results may not be favorable when the sensor of interest is not (or not concurrent with) the one from which the profiles are derived [26,32]. Overcoming this temporal limitation, profiles from global reanalysis data provide a flexible temporal resolution (typically 3 and 6 hourly) and are a practical alternative to the radiosonde’s spatial constraint [21]. Reanalysis datasets are global gridded extended homogeneous time series, with no spatial or temporal gaps [41].

Barsi et al. (2005, 2003) [27,29] proposed an atmospheric web-based correction tool (Atmospheric Correction Parameter Calculator—ACPC) for Landsat 5, 7, and 8 thermal bands. It uses reanalysis profiles from the National Centers for Environmental Prediction (NCEP) and the code of the moderate resolution atmospheric transmission (MODTRAN) model [42] to directly provide the three atmospheric parameters for AC. Additionally, researchers have evaluated and compared the efficacy of different reanalysis and satellite-based profile products for AC/LST retrieval. NCEP reanalyses (Reanalysis 1 [43] and FNL [44]) and a moderate resolution radiometer (MODIS) atmospheric profiles product (MOD07) [45] were analyzed for different sites in Spain [24,33,46] and in China [47], for TIR sensors such as Landsat, ASTER, and HJ-1B IRS. Coll et al. (2012) [32], in Spain, added the satellite-based profiles from Atmospheric Infrared Sounder (AIRS) [48] to the comparison. Rosas et al. (2017) [21] also included the European Centre for Medium-Range Weather Forecast’s (ECMWF) ERA-Interim reanalysis product [49] along with previous profiles for Landsat 8 LST estimates in Saudi Arabia. More recently, assessments with a greater number of different profile products were carried out. Meng and Cheng (2018) [20] evaluated reanalysis products from Modern-Era Retrospective Analysis for Research and Applications (MERRA) [50] 1 and 2, ERA-Interim, NCEP FNL and NCEP Reanalysis 2 [51], and Japanese Reanalysis (JRA-55) [52]. Their validation included profiles distributed around the globe and 32 Landsat 8 band 10 images over China. Yang et al. (2020) [40] assessed seven profile sources (AIRS, MOD07, ERA-Interim, MERRA 2, and NCEP’s FNL, Reanalysis 2, and GFS) using 17 radiosonde stations over Europe and in situ LST measurements in China. Overall, the accuracies of AC parameters and retrieved LST data using satellite-derived profiles were lower than those using reanalysis methods [32,33,40,46,47].

However, reanalysis profiles also have their disadvantages. The spatial resolution of several degrees (varying for each product) can be considered low. The accuracy is usually poorer for regions with less coverage of permanent observatories, such as the oceans and many Southern Hemisphere zones [41,53,54]. Since the data are spaced at grid points with time intervals commonly of 6 h, the description of meteorological phenomena on a sub-grid or variable time scale may be affected [55]. On the other hand, modern numerical weather models benefit from computing performance and physical processes parameterization in downscaling the reanalysis data. Mesoscale atmospheric models use global (re)analysis models as initial and boundary conditions for local applications [56,57]. Lee et al., (2020) [58] employed high-resolution (1.5 and 12 km) Numerical Weather Prediction (NWP) models distributed by the Korea Meteorological Administration (KMA) as input atmospheric data for AC and sea surface temperature estimation with visible infrared imaging radiometer suite (VIIRS) bands. Their KMA NWP dataset is restricted to eastern Asia, but the approach of using high-resolution NWP models combined with RTM provides an interesting background for studies of TIR remote sensing atmospheric correction.

The Weather Research and Forecasting (WRF) model [59] is an atmospheric modeling system designed for both research and NWP. It is a state-of-the-art mesoscale model and the world’s most widely used model for these purposes. Non-hydrostatic, open-source, free, community-based, and with a wide range of parameterization options, the WRF model provides a spectrum of capabilities for a variety of applications in atmospheric science and weather prediction [60]. Hence, the WRF model has been extensively employed for estimating high-resolution meteorological data [61,62,63,64,65,66].

It is imperative to assess whether the use of the WRF model to generate high-resolution atmospheric vertical profiles, in conjunction with an RTM, results in a higher AC/LST retrieval accuracy. Moreover, almost all studies that evaluated different profile sources for LST retrieval were performed over Asia and Europe. To the best of our knowledge, no study has carried out such an assessment in South America. There is also a lack of studies using newer and better reanalysis profiles (e.g., ERA5 [67] and NCEP Climate Forecast System Version 2 (CFSv2) [68]) for AC and LST estimation [38,69].

This study conducted simulations with the WRF model using NCEP CFSv2 reanalysis data as initial and boundary conditions. Its objective was to generate high-resolution vertical profiles, improving the spatial, temporal, and vertical resolutions of the global reanalysis. The intention was to investigate the utility of these profiles in TIR atmospheric correction and LST retrieval, in relation to the ACPC web-tool and profiles extracted directly from the NCEP CFSv2 reanalysis. We used Landsat 8 TIRS band 10 as an example to retrieve LST values through an RTE inversion-based algorithm in conjunction with the MODTRAN radiative transfer model. The accuracy assessment was performed using local radiosonde observations in Southern Brazil.

2. Materials and Methods

2.1. Study Area and In Situ Radiosonde Data

The Porto Alegre International Airport (SBPA), Rio Grande do Sul State, Brazil was selected as the study area. The airport includes a radiosonde station, which made this site a useful environment for studies that aimed to evaluate atmospheric profiles. The selected area covers the official limits of the Anchieta district [70], with an area of around 9.2 km² (Figure 1). According to Köppen’s climate classification, the local climate is subtropical humid with hot summers (Cfa) [71], with a mean annual air temperature of 19.6 °C, mean annual relative humidity of 76.1%, and mean annual precipitation of 1397 mm [72].

The SBPA station is located at 30.00° S and 51.18° W, 3.0 m above mean sea level, and its identifier number is 83971. In this station, radiosondes are launched twice a day, at 00:00 and 12:00 UTC. In this study, we used the 12:00 UTC radiosonde profiles, as this is the closest time to the Landsat 8 crossing time over the study site (~13 UTC). Data were obtained from the University of Wyoming website (http://weather.uwyo.edu/upperair/sounding.html, accessed on 28 October 2021). This dataset allows the vertical structure of the atmosphere to be characterized, with profiles of air temperature, pressure, and humidity with up to 99 vertical levels. The radiosonde observations, as well as the parameters calculated from them, are considered as ground truth for the assessments in this study.

2.2. Landsat 8 Satellite Data and Case Days

The Landsat mission has been providing moderate-resolution space-based surface observations for almost 50 years. Landsat 8 is the most recent operational satellite of the series and was launched in February 2013. It carries a two-sensor payload: the Operational Land Imager (OLI), which has nine reflective (visible, near-infrared, and short-wave infrared) bands with a 30 m spatial resolution and the Thermal Infrared Sensor (TIRS) with two bands in the TIR region. The TIRS bands have a 100 m spatial resolution but it is resampled and provided at 30 m by the United States Geological Survey (USGS), to be consistent with the OLI bands [73,74].

In this study, we acquired all Landsat 8 images (Collection 1) available under daily clear-sky conditions over the study area from 2013 to 2019. This resulted in a total of 27 scenes at Path-Row 221-81 (Figure 1), with an acquisition scene center time around 13:20 UTC. The full-swath Landsat data were reduced to a subset for a 10,184-pixel region covering the study area of Figure 1.

Table 1. Information about Landsat 8 data used in the paper (27 images from 2013 to 2019, Path-Row 221-81). The last columns present meteorological conditions measured by the SBPA station.

Scene ID	Acquisition Date	Case Day	Season	Air T [°C]	PWV [g/cm2]
LC82210812013322LGN01	18 November 2013	1	spring	23	1.57
LC82210812013338LGN01	4 December 2013	2	spring	21	2.43
LC82210812014037LGN01	6 February 2014	3	summer	31	4.10
LC82210812014293LGN01	20 October 2014	4	spring	19.8	1.18
LC82210812014341LGN01	7 December 2014	5	spring	26.6	1.99
LC82210812015024LGN01	24 January 2015	6	summer	24.4	2.66
LC82210812015056LGN01	25 February 2015	7	summer	24.6	3.26
LC82210812015312LGN01	8 November 2015	8	spring	24.4	2.79
LC82210812016075LGN01	15 March 2016	9	summer	22.8	2.23
LC82210812016235LGN02	22 August 2016	10	winter	9.4	0.91
LC82210812016347LGN01	12 December 2016	11	summer	23.4	3.27
LC82210812017093LGN01	3 April 2017	12	autumn	22	3.31
LC82210812017173LGN00	22 June 2017	13	winter	10.6	1.85
LC82210812017205LGN00	24 July 2017	14	winter	13.4	1.93
LC82210812017237LGN00	25 August 2017	15	winter	22.2	2.82
LC82210812017317LGN00	13 November 2017	16	spring	20	1.83
LC82210812017349LGN00	15 December 2017	17	spring	25.2	3.17
LC82210812018048LGN00	17 February 2018	18	summer	23.4	2.36
LC82210812018112LGN00	22 April 2018	19	autumn	19	3.26
LC82210812018160LGN00	9 June 2018	20	autumn	6.6	0.84
LC82210812018240LGN00	28 August 2018	21	winter	9.4	0.60
LC82210812018272LGN00	29 September 2018	22	spring	23	3.60
LC82210812018320LGN00	16 November 2018	23	spring	23.6	1.74
LC82210812019083LGN00	24 March 2019	24	autumn	25.8	2.53
LC82210812019099LGN00	9 April 2019	25	autumn	19.6	1.72
LC82210812019227LGN00	15 August 2019	26	winter	10.8	1.08
LC82210812019323LGN00	19 November 2019	27	spring	24.8	2.33

presents the specifications of the Landsat 8 data utilized in the study, as well as the case days that henceforward will be used to refer to each dataset. Seeking to illustrate the meteorological conditions near the acquisition time, we also include in Table 1 the air temperature and the water vapor content in the atmosphere, using the precipitable water vapor (PWV) data for the entire radiosonde column measured by the SBPA station. TIRS band 10 (10.60–11.19 µm) was used for RTE-based LST retrieval (Section 2.6.2) and OLI bands 4—red (0.64–0.67 µm) and 5—near-infrared (0.85–0.88 µm) for land surface emissivity estimation (Section 2.5).

2.3. Reanalysis Data

The NCEP Climate Forecast System version 2 (CFSv2) [68] reanalysis data are produced using the NCEP Global Forecasting System (GFS) atmospheric model and the Gridpoint Statistical Interpolation (GSI) analysis system with three-dimensional variational data assimilation (3D-Var). This is arranged in grids with a horizontal resolution of 0.5° × 0.5° and in 37 vertical (pressure) levels (1000—1 mbar), with 0.205° for surface parameters. We used CFSv2 reanalysis data from the 6-hourly product as initial and boundary conditions for the WRF simulations. In addition, profiles retrieved directly from NCEP CFSv2 were included in the analysis, to assess the WRF model downscaling performance. These profiles were extracted from the grid point closest to the SBPA station.

2.4. WRF Model Configuration

The WRF model version 4.1.2 with the advanced research WRF (ARW) dynamical solver [59,75] was used to perform high-resolution numerical simulations. We configured the WRF domains with two nested grids, in one-way mode, centered at the SBPA station, with horizontal resolutions of 12 km (G12) and 3 km (G03) (4:1 parent grid ratio) (see Figure 1) and 33 sigma vertical levels with 50 hPa as the top pressure value. The time step used for the outermost domain was 72 s with a 4:1 parent ratio. Geographical static data such as land use, topography, albedo, and reflectance, were inserted into the modeling process using the land use categories of USGS.

Some atmospheric physical processes cannot be directly resolved by the numerical model and so need to be parameterized. Parameterization schemes represent the contribution of unresolved but important phenomena in terms of variables resolved at the model discrete grid [76,77,78]. The WRF model has a vast range of parameterization scheme options. The physics parameterization chosen for our simulations included: Purdue Lin microphysics parameterization [79]; the Yonsei University (YSU) Planetary Boundary Layer (PBL) scheme [80]; the Betts–Miller–Janjic (BMJ) cumulus scheme [81]; Dudhia shortwave radiation [82]; the Rapid Radiative Transfer Model (RRTM) longwave radiation option [83]; the Unified NOAH Land-Surface Model (LSM) [84]; the revised MM5 surface-layer scheme [85]. In Diaz et al. (2021) [86], some parametrization options were tested for the same SBPA area, and the set used here is based on these results and those of Santos and Nascimento (2016) [87]. The WRF configurations are summarized in

Table 2. Overview of WRF model setting.

WRF Model Configuration
Version	4.1.2
Dynamical solver	ARW
Boundary conditions	NCEP CFSv2
Map projection	Lambert
Grid size	Domain 1: (119 × 116) × 33 Domain 2: (169 × 165) × 33
Horizontal resolution	Domain 1: 12 km Domain 2: 3 km
Nesting	One-way
Time step	72 s
Static geographical data	USGS
Cloud Microphysics	Purdue Lin
Planetary Boundary Layer (PBL)	Yonsei University (YSU)
Cumulus	Betts–Miller–Janjic (BMJ) 1
Shortwave Radiation	Dudhia
Longwave Radiation	Rapid Radiative Transfer Model (RRTM)
Land Surface Model (LSM)	Unified NOAH
Surface-layer	Revised MM5

¹ Domain 1 only.

[88].

The WRF model was run using the above configurations for each of the case days in Table 1. We conducted simulations of 24 h duration starting at 00:00 UTC and the resulting profiles were extracted at 12:00 UTC for the grid point closest to the SBPA station, to match with the local radiosonde observations. Hence, the first 12 h of the simulation was considered to be spin-up time. The model output for each domain was stored every 30 min.

2.5. Land Surface Emissivity Estimation

Land surface emissivity (LSE) is one of the key parameters for retrieving LST from remote sensing data. It is a measure of the surface’s intrinsic ability to convert heat energy into radiant energy. LSE is a function of the composition, roughness, and moisture of the surface and the observation conditions [74,89,90]. Among the methods for LSE retrieval from space, those based on the normalized difference vegetation index (NDVI) are operational and are the most frequently applied, with satisfactory results [89,91,92,93,94,95,96]. Sekertekin and Bonafoni (2020a, 2020b) [74,97] examined the influence of six NDVI-based LSE models on the performance of LST retrieval. Based on their results, we opted to use the NDVI threshold method (NDVI^THM) of Sobrino et al. (2008) [91]. The LSE estimation from satellites inevitably has errors and NDVI-based methods also have their limitations, especially in urban environments [98]. However, a recent study [99] reported that the use of LSE data estimated by different methods resulted in no significant variations in the LST accuracy. Nevertheless, assessing the LSE estimation is beyond the scope of this paper.

To calculate the NDVI from Landsat 8 data, the first step is to convert the digital number (DN) values of bands 4 (red) and 5 (near-infrared—NIR) to reflectances using Equation (1) [97,100]:

$${\mathsf{\rho }}_{\mathsf{\lambda }}=\frac{{\mathrm{M}}_{\mathrm{p}}·{\mathrm{Q}}_{\mathrm{CAL}}+{ \mathrm{A}}_{\mathrm{p}}}{\sin {\mathsf{\theta }}_{\mathrm{SE}}}$$

(1)

where ${\mathsf{\rho }}_{\mathsf{\lambda }}$ is the reflectance of the corresponding band, M_P is the multiplicative rescaling factor of the corresponding band, Q_CAL is the calibrated and quantized standard product of pixel values (DNs), A_P is the additive rescaling factor of the corresponding band, and ${\mathsf{\theta }}_{\mathrm{SE}}$ is the local sun elevation angle.

Then, the NDVI is obtained from Equation (2):

$$\mathrm{NDVI}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{ \mathsf{\rho }}_{\mathrm{R}}}{{\mathsf{\rho }}_{\mathrm{NIR}}+{ \mathsf{\rho }}_{\mathrm{R}}}$$

(2)

where ${\mathsf{\rho }}_{\mathrm{NIR}}$ is the reflectance of the NIR band and ${\mathsf{\rho }}_{\mathrm{R}}$ is the reflectance of the red band. It is not necessary to correct the atmospheric effects in the red and NIR bands to estimate LSE [101].

From the NDVI, it is possible to calculate the fractional vegetation cover (${\mathrm{P}}_{\mathrm{V}}$) from Equation (3) [102]:

$${\mathrm{P}}_{\mathrm{V}}={\left[\frac{\mathrm{NDVI}-{ \mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{ \mathrm{NDVI}}_{\min }}\right]}^2$$

(3)

where ${\mathrm{NDVI}}_{\min }$ = 0.2 and ${\mathrm{NDVI}}_{\max }$ = 0.5 in a global context [30,91,103]. The ${\mathrm{P}}_{\mathrm{V}}$ is an important factor in the LSE estimation.

The NDVI^THM proposed by (Sobrino et al., 2008) [91] estimates the LSE considering three different cases as presented in Equation (4), for Landsat 8 [74]:

$$\mathsf{\epsilon }=\{\begin{array}{c}0.979-0.035{\mathsf{\rho }}_{\mathrm{R}}\\ 0.004{\mathrm{P}}_{\mathrm{V}}+0.986\\ 0.99\end{array}\begin{array}{ccc}\hfill \hspace{1em}& \hfill \mathrm{NDVI}& \hfill < 0.2\\ \hfill \hspace{1em}0.2 \le & \hfill \mathrm{NDVI}& \hfill \le 0.5\\ \hfill \hspace{1em}& \hfill \mathrm{NDVI}& \hfill > 0.5\end{array}$$

(4)

where $\mathsf{\epsilon }$ is the land surface emissivity (LSE). For NDVI < 0.2, the pixel is considered to be 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 to be composed of a mixture of bare soil and vegetation and the LSE depends on the ${\mathrm{P}}_{\mathrm{V}}$ value. The pixels with NDVI values higher than 0.5 are considered to be fully vegetated areas and the emissivity is assumed to be 0.99.

2.6. Atmospheric Correction and LST Retrieval

2.6.1. Atmospheric Parameters Calculation with MODTRAN and ACPC

MODTRAN (moderate resolution atmospheric transmission) is a commercial and widely employed atmospheric RTM developed by the U.S. Air Force and Spectral Sciences Inc. [42]. The present study used the MODTRAN4 v3r1 [104] 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; (iv) WRF G03. The atmospheric parameters calculated from SBPA profiles were treated as ground truth data.

The methodology of Barsi et al. (2003) [29] was adopted to fulfill the profiles. To predict space-reaching atmospheric parameters, MODTRAN requires atmospheric profiles reaching “space”, or 100 km above sea level. Since the radiosondes and NCEP CFSv2 stretch from the surface to about 30 and 50 km, respectively, the upper atmosphere layers (to 100 km) were extracted from the MODTRAN standard atmospheres and pasted into our site-specific profiles. We used the standard mid-latitude summer profile [105] for case days in hot seasons (spring and summer) and the mid-latitude winter profile [106] for those in cold seasons (autumn and winter). The WRF-simulated profiles refined the NCEP CFSv2 profiles by increasing the number of levels in the portion of the atmosphere closest to the surface. The model was set up with the top pressure level at 50 hPa (~18 km). Thus, these profiles were completed first with the remaining levels from the reanalysis and then, starting at 50 km, with MODTRAN standard atmospheres.

This approach results in surface-to-space vertical profiles of air temperature, pressure, and water vapor. The profiles constructed for MODTRAN from the NCEP CFSv2 reanalysis count 40 (43) vertical levels in the warm (cold) seasons. Using WRF, the number of vertical levels in the filled profiles increases to 45 (48). The radiosondes describe the atmospheric profile in more detail, but the number of vertical levels varies, averaging 88 and reaching 109 in our final profiles. These completed profiles were inserted into a MODTRAN input file and then processed [27,29].

MODTRAN outputs are provided in the model’s spectral resolution, so an integration must be performed between the bounds defined by the spectral response curve of the sensor band (Equation (5)) [4]:

$$\mathrm{var} \left({\mathsf{\lambda }}_{\mathrm{i}}\right)=\frac{{\int }_{{\mathsf{\lambda }}_{\mathrm{i},\min }}^{{\mathsf{\lambda }}_{\mathrm{i},\max }}\mathrm{var} \left(\mathsf{\lambda }\right){\mathrm{R}}_{\mathrm{s}}\left(\mathsf{\lambda }\right)\mathrm{d}\mathsf{\lambda }}{{\int }_{{\mathsf{\lambda }}_{\mathrm{i},\min }}^{{\mathsf{\lambda }}_{\mathrm{i},\max }}{\mathrm{R}}_{\mathrm{s}}\left(\mathsf{\lambda }\right)\mathrm{d}\mathsf{\lambda }}$$

(5)

where ${\mathrm{R}}_{\mathrm{s}}$ is the spectral response of the sensor at the center wavelength ${\mathsf{\lambda }}_{\mathrm{i}}$ of a band with a spectral window of ${\mathsf{\lambda }}_{\mathrm{i},\min }$–${\mathsf{\lambda }}_{\mathrm{i},\max }$ and $\mathrm{var}$ alternatively denotes the atmospheric parameter (transmittance, upwelling, or downwelling radiance) value extracted from the MODTRAN output file, between the wavelengths ${\mathsf{\lambda }}_{\mathrm{i},\min }$ and ${\mathsf{\lambda }}_{\mathrm{i},\max }$. In this case, the spectral response curve is from the Landsat 8 TIRS band 10 (TIRS10).

Additionally, we included in the comparative analysis the atmospheric parameters estimated by NASA’s well-established Atmospheric Correction Parameter Calculator (ACPC) web-tool [27,29]. As mentioned above, the ACPC uses NCEP reanalysis profiles (with 1° × 1° horizontal resolution and 28 vertical levels), MODTRAN code, and a suite of integration algorithms to provide the AC parameters for the particular date, time, and location inputted. The reanalysis profiles used by ACPC are to about 30 km, so they are fulfilled using MODTRAN standard atmospheres. Completed ACPC profiles have 33 vertical levels. The mid-latitude standard upper profiles varied according to the season of each case day. The option of using the atmospheric profile from the closest integer coordinate to the inputted location (SBPA station) was set.

2.6.2. Radiative Transfer Equation (RTE)-Based LST Retrieval Method

The inverse solution of the RTE [28] is a direct and a priori the most appropriate procedure for LST retrieval using a single TIR band [101]. The RTE applied to a particular TIR band/wavelength ($\mathsf{\lambda }$) can be simplified and given by:

$${\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}=\left[{\mathsf{\epsilon }}_{\mathsf{\lambda }}{\mathrm{B}}_{\mathsf{\lambda }}\left({\mathrm{T}}_{\mathrm{s}}\right)+\left(1-{\mathsf{\epsilon }}_{\mathsf{\lambda }}\right){\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }\right]{\mathsf{\tau }}_{\mathsf{\lambda }}+{\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }$$

(6)

where ${\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}$ (W·m⁻²·sr⁻¹·µm⁻¹) is the at-sensor (TOA) spectral radiance of the corresponding TIR band (in this paper, TIRS10), ${\mathsf{\epsilon }}_{\mathsf{\lambda }}$ refers to the LSE (dimensionless), ${\mathrm{B}}_{\mathsf{\lambda }}$ (W·m⁻²·sr⁻¹·µm⁻¹) is the black body radiance, ${\mathrm{T}}_{\mathrm{s}}$ (Kelvin) represents the LST, ${\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }$ and ${\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }$ (W·m⁻²·sr⁻¹·µm⁻¹) refer to the downwelling and upwelling radiances, respectively, and ${\mathsf{\tau }}_{\mathsf{\lambda }}$ is the atmospheric transmittance (dimensionless). Hence, the emitted radiance for a black body at a temperature ${\mathrm{T}}_{\mathrm{s}}$ is given by the inversion of Equation (6):

$${\mathrm{B}}_{\mathsf{\lambda }}\left({\mathrm{T}}_{\mathrm{s}}\right)=\frac{{\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}-{\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }-\mathsf{\tau }\left(1-{\mathsf{\epsilon }}_{\mathsf{\lambda }}\right){\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }}{{\mathsf{\tau }}_{\mathsf{\lambda }}{\mathsf{\epsilon }}_{\mathsf{\lambda }}}$$

(7)

where the ${\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}$ of Landsat 8 TIRS10 is obtained by converting the DN values (${\mathrm{Q}}_{\mathrm{CAL}}$) applying Equation (8) [100]:

$${\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}={ \mathrm{M}}_{\mathrm{L}}·{\mathrm{Q}}_{\mathrm{CAL}}+{ \mathrm{A}}_{\mathrm{L}}$$

(8)

where ${\mathrm{M}}_{\mathrm{L}}$ and ${\mathrm{A}}_{\mathrm{L}}$ are the multiplicative and additive rescaling factors of the corresponding band (TIRS10). In addition, ${\mathrm{B}}_{\mathsf{\lambda }}$ comes from Planck’s law:

$${\mathrm{B}}_{\mathsf{\lambda }}\left(\mathrm{T}\right)=\frac{{\mathrm{C}}_1}{{\mathsf{\lambda }}^5\left({\mathrm{e}}^{{\mathrm{C}}_2/\mathsf{\lambda }\mathrm{T}}-1\right)}$$

(9)

where ${\mathrm{C}}_1$ = 1.19104 × 10⁸ W·μm⁴·m⁻²·sr⁻¹ and ${\mathrm{C}}_2$ = 14,387.7 µm·K are Planck’s radiation constants. Thus, ${\mathrm{T}}_{\mathrm{s}}$ is calculated by inverting Planck’s law in Equation (7):

$${\mathrm{T}}_{\mathrm{s}}=\frac{{\mathrm{C}}_2}{\mathsf{\lambda }}{\left[\ln \left(\frac{{\mathrm{C}}_1}{{\mathsf{\lambda }}^5\left[\frac{{\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}-{\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }}{{\mathsf{\epsilon }}_{\mathsf{\lambda }}{\mathsf{\tau }}_{\mathsf{\lambda }}}-\left(\frac{1-{\mathsf{\epsilon }}_{\mathsf{\lambda }}}{{\mathsf{\epsilon }}_{\mathsf{\lambda }}}\right){\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }\right]}+1\right)\right]}^{-1}$$

(10)

Finally, Equation (10) can be simplified and the LST (${\mathrm{T}}_{\mathrm{s}}$) from Landsat 8 TIRS10 is estimated as:

$$\mathrm{LST}=\frac{{\mathrm{K}}_2}{\ln \left(\frac{{\mathrm{K}}_1}{\frac{{\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}-{\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }-{\mathsf{\tau }}_{\mathsf{\lambda }}\left(1-{\mathsf{\epsilon }}_{\mathsf{\lambda }}\right){\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }}{{\mathsf{\tau }}_{\mathsf{\lambda }}{\mathsf{\epsilon }}_{\mathsf{\lambda }}}}+1\right)}$$

(11)

where ${\mathrm{K}}_1$ and ${\mathrm{K}}_2$ refer to calibration constants, whose values for Landsat 8 TIRS10 are 774.89 W·m⁻²·sr⁻¹·µm⁻¹ and 1321.08 K, respectively [100]. Henceforth, the spectral notation ($\mathsf{\lambda }$) will be omitted, since here only a single TIR band is used.

The aforementioned procedure was applied with τ, ${\mathrm{L}}^{\uparrow }$, and ${\mathrm{L}}^{\downarrow }$ calculated using profiles from:

• SBPA local radiosonde;
• NCEP CFSv2 reanalysis;
• WRF G12;
• WRF G03;
• ACPC.

The LST images retrieved using the atmospheric parameters from SBPA radiosondes are considered as “ground truth”, and LST retrievals considering the other cases were compared to this, as detailed in the next section.

2.7. Metrics for Performance Evaluation

To evaluate the performance of the WRF model and other profiles we take into account the atmospheric parameters (τ, ${\mathrm{L}}^{\uparrow }$, and ${\mathrm{L}}^{\downarrow }$) 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, the 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 [31,107,108]. The metrics’ equations are given as follows:

$$\mathrm{R}=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{p}}_{\mathrm{i}}- \overline{\mathrm{p}}\right)\left({\mathrm{o}}_{\mathrm{i}}- \overline{\mathrm{o}}\right)}{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{p}}_{\mathrm{i}}- \overline{\mathrm{p}}\right)}^2}\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{o}}_{\mathrm{i}}- \overline{\mathrm{o}}\right)}^2}}$$

(12)

$$\mathrm{bias}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{p}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right)$$

(13)

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{p}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right|$$

(14)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{p}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right)}^2}$$

(15)

where ${\mathrm{p}}_{\mathrm{i}}$ and ${\mathrm{o}}_{\mathrm{i}}$ are pairwise predicted and observed values, respectively, $\mathrm{n}$ is the number of pairs, and terms with overbars are respective mean values.

R scores range from −1 to 1, and values approaching 1 indicate stronger correlations. Bias is useful in determining if the model is underestimating (bias < 0) or overestimating (bias > 0) the observed values. MAE is the average of the absolute difference between the predictions and the observations, despite model overestimation or underestimation. It is a negatively oriented index (values closer to zero are better). RMSE is also negatively oriented. It indicates the deviation between modeled and observed values and is more sensitive to outliers as the errors are squared before summing.

3. Results and Discussion

3.1. Evaluation of Atmospheric Parameters

3.1.1. Overall Results

The AC parameters (τ, ${\mathrm{L}}^{\uparrow }$, and ${\mathrm{L}}^{\downarrow }$) calculated with different sources of estimated vertical profiles (CFSv2, WRF G12, WRF G03, and ACPC) were compared against those using observational SBPA radiosondes (Figure 2). In

Table 3. Statistical metrics of atmospheric parameters estimated from different atmospheric profiles. Bias, MAE, and RMSE are in the units of the parameter of interest, i.e., W·m⁻²·sr⁻¹·µm⁻¹ for atmospheric radiances and the transmittance is dimensionless.

		CFSv2	WRF G12	WRF G03	ACPC
Transmittance	R	0.96	0.93	0.93	0.97
	bias	0.01	0.00	0.00	0.01
	MAE	0.02	0.03	0.03	0.02
	RMSE	0.03	0.04	0.04	0.03
Upwelling	R	0.97	0.94	0.94	0.98
	bias	−0.12	−0.04	−0.02	−0.06
	MAE	0.21	0.24	0.24	0.16
	RMSE	0.27	0.33	0.34	0.20
Downwelling	R	0.97	0.95	0.95	0.98
	bias	−0.15	−0.05	−0.03	0.08
	MAE	0.28	0.30	0.30	0.21
	RMSE	0.35	0.42	0.43	0.27

, the accuracies of atmospheric parameter estimations are presented. All the profile sources provide AC parameters with high correlation coefficients, all greater than 0.9 in relation to the reference (SBPA). The R values of ACPC are slightly better, followed by those of CFSv2. There is a general but small tendency to overestimate 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 were from CFSv2. With respect to MAE and RMSE, the best results were from ACPC followed by CFSv2. The largest RMSE value found was 0.43 W·m⁻²·sr⁻¹·µm⁻¹ for the downwelling calculated with WRF G03 profiles.

The higher negative bias of the atmospheric radiances with NCEP CFSv2 in our findings may be due to the fact that these reanalysis profiles have the lowest level at 1000 hPa, which corresponds to around 60–250 m for the SBPA station in our case days. Therefore, the lowest layer of the atmosphere (which typically presents the largest water vapor content and warmest temperature) is neglected in these profiles. This fact was indicated by Coll et al. (2012) [32] and Meng and Cheng (2018) [20]. We tried to reduce this limitation by downscaling the CFSv2 reanalysis profiles with the WRF model. The WRF profiles bring the first level to around 1 m above the surface.

Figure 3 illustrates the performance of the AC parameters from the different profile sources via Taylor diagrams [109]. These quantify the agreement of the modeled profiles with the observational profiles in terms of graphical representations of the three statistics combined: standard deviation, root mean squared deviation (RMSD, the same as RMSE), and correlation coefficient. Figure 3 corroborates the results shown in Table 3. It indicates that ACPC obtained better results than the other profiles for all three atmospheric parameters. Next was the CFSv2 reanalysis profiles. The diagrams clearly point out that the profiles from WRF G12 and G03 had very similar accuracies. In fact, no significant statistical differences were found between the parameters from the WRF grids G12 and G03. This 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 [57,110,111,112].

Our RMSE range is in agreement with the findings for reanalysis profiles in Meng and Cheng (2018) [20] 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) [40] were overall slightly lower than ours. However, their assessment only includes profiles from Europe. In Diaz et al. (2021) [85], the vertical distributions of air temperature and water vapor for the SBPA radiosondes, CFSv2 reanalysis, and WRF (G12 and G03) profiles were compared. The evaluation of the profiles themselves showed that both CFSv2 and WRF models skillfully represent the vertical profiles of temperature and water vapor. Nevertheless, the statistical metrics indicated that increasing the horizontal resolution did not significantly improve the quality of the simulated atmospheric profile. This is in line with the above results for the atmospheric parameters comparison.

3.1.2. Analysis by Meteorological Conditions

The errors in the estimation of atmospheric parameters for each case day are shown in Figure 4. This shows that, despite the fact that ACPC presented the best overall metrics, none of the profile sources performed best in all cases. For instance, in eight of the case days, one of the WRF profiles had the best results in calculating the downwelling radiance. For case day 23, the WRF model plainly produced the largest errors, whereas for case days such as day 26, the model successfully reduced the largest error in the driving reanalysis data. Therefore, it is pertinent to assess the dependence of the performances of different profile sources on the meteorological conditions.

Figure 5 shows the RMSE values of the AC parameters calculated taking into account specific ranges of precipitable water vapor (PWV) for the entire atmospheric column (measured by SBPA radiosondes) and air temperature observed at the SBPA station (Table 1). The WRF model performs better under dry and cold conditions. It outperforms the other profile sources in the range of 2–3 g/cm² PWV. Under almost all moisture and temperature conditions, the ACPC profiles presented better results than the NCEP CFSv2 profiles. All the different profiles presented minor errors in situations of the lowest water vapor and air temperature. However, it is important to mention that the number of case days under these conditions was less than for the other ranges. This is consistent with the research of Meng and Cheng (2018) [20], where it was noted that the best global RMSEs were for profiles with a water vapor content lower than 1 g/cm². Furthermore, they indicated that the largest RMSEs occurred when the water vapor was between 3 and 4 g/cm². With respect to the air temperature, our results suggest that errors in estimating the atmospheric parameters are higher on warmer days.

3.2. Application to RTE-Based LST Retrieval

To further assess the different atmospheric profiles, the LSTs retrieved by RTE inversion with atmospheric parameters, Landsat 8 TIRS10 radiance, and NDVI^THM emissivity were compared with each other. 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 (11) were the same for each pixel of the scenes. Therefore, the differences between LST values that were retrieved with SBPA profiles and those from other sources were due to the discrepancies among the profiles [40]. Figure 6 illustrates an example (case day 25) of the LST spatial distribution in the study area. Note that the hottest pixels are over the airstrip and adjacent buildings of the airport.

Histograms of LST errors for all the 10,184 pixels in the study area of the 27 case days are shown in Figure 7. 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 in the range of ±1 K. Yang et al. (2020) [40] also found the most LST differences in this range, using different reanalysis and satellite-derived profiles. The histograms in Figure 7b,c indicate that WRF profiles tend to overestimate the LST, whereas ACPC tends to underestimate it (Figure 7d). Using WRF profiles, LST errors can reach to more than 4 K, although only 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 4. Statistical metrics of land surface temperatures (LST) retrieved using atmospheric parameters from different profile sources. The LST values calculated with SBPA parameters were considered as reference values.

		CFSv2	WRF G12	WRF G03	ACPC
LST [K]	R	0.99	0.99	0.99	0.99
	bias	0.23	0.32	0.36	−0.38
	MAE	0.54	0.79	0.81	0.56
	RMSE	0.55	0.79	0.82	0.56

summarizes the metrics of the LST retrieval comparative analysis. Overall, the LST values of the four profile sources analyzed in this study were found to be 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 bias. This corresponds with the histogram analysis in Figure 7. The mean error criteria (MAE and RMSE) indicate that the profiles with the best performance in the AC- and RTE-based LST retrieval are, in descending order: CFSv2, ACPC, WRF G12, and WRF G03. The differences between CFSv2 and ACPC overall MAE and RMSE values were very small, and similarly for WRF G12 and G03.

Comparing these results with previous studies that evaluated the application of different atmospheric profiles for LST retrieval, Meng and Cheng (2018) [20] reported overall LST RMSE values higher than ours for eight different reanalysis profile sources analyzed around the globe. All their eight average RMSEs were larger than 1 K. The authors mention that ERA-Interim and MERRA (6 h product) showed the lowest RMSEs, at 1.09 K and 1.07 K respectively. They also indicated an average tendency to overestimate the LST, except for JRA-55. In Yang et al. (2020) [39], RMSEs were smaller than 0.6 K using profiles from ERA-Interim, MERRA2, NCEP GFS, and NCEP FNL, over Europe. Their worst accuracy (RMSE of 1 K) was obtained using MOD07 satellite-based profiles. Retrieving LST from three MODIS thermal bands, Pérez-Planells et al. (2015) [33] 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 8 displays the distribution of LST bias and RMSE among the case days. Figure 8a evidences the ACPC and WRF settings’ average tendency to 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. However, the metrics of mean errors are greater for WRF profiles, as in 11 case days the model improves the results of the reanalysis. Furthermore, the finer grid (G03) succeeded in downscaling the G12 11 times. In fact, the largest errors were achieved when using the WRF model (e.g., for case days 23, 19, and 17), contributing to the higher overall RMSE values. In general, the days with larger errors in atmospheric parameters (Figure 4) were the days with larger LST RMSEs, as in case day 23. Conversely, for case day 22, the errors in atmospheric parameters using ACPC were less than those using WRF profiles, but the highest LST RMSE on this day was with ACPC. Jiménez-Muñoz et al. (2010) [24] suggest that cases like this may be explained by compensation among the AC parameter errors, for instance, a significant positive difference in transmittance and significant but negative differences in the atmospheric radiances. Figure 9 shows how errors in the atmospheric parameters propagate to the retrieved LST. The sensitivity analysis of Sekertekin and Bonafoni (2020) [74] reported that an uncertainty of ±0.01 in atmospheric transmittance had an impact of ±0.97 K on RTE-based LST retrieval. For upwelling and downwelling, uncertainties of ±10% led to ±1.82 K and ±0.07 K errors in LST, respectively. Rosas et al. (2017) [21] mentioned that an introduced uncertainty of 20% in the relative humidity profile could result in LST errors as large as 1.5 K in arid environments.

In summary, the attempt at 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. Using WRF simulations, we improved the vertical resolution in the lowest atmospheric levels. However, no significant improvement was found in the AC using the WRF profiles. In some cases, using a finer grid resolution profile resulted in even greater differences in atmospheric parameters and LST estimation. Rosas et al. (2017) [21] 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 data with higher resolution naturally tend to represent the atmospheric parameters better, it is not a strictly direct relationship [20,40]. Furthermore, the ACPC, which uses profiles with 1° × 1° (~111 km) horizontal resolution, showed good results. Although these profiles have a coarser horizontal resolution, previous studies have found satisfactory results using the ACPC, even surpassing other methods [24,31,46,93,113,114,115,116,117]. It is important to note that in this study the WRF profiles were extracted at 12 UTC to match the available radiosonde data. Nevertheless, the extraction could be set for the exact time of the satellite overpass. In the ACPC, and in general for reanalysis profiles, this time synchronization is achieved through linear interpolation, which may not be the most appropriate strategy for sampling weather fronts and diurnal heating cycles [20,27,29].

4. 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 profiles were incorporated into the MODTRAN RTM to compute the atmospheric correction parameters, which were then applied in the RTE to retrieve LST values from Landsat 8 TIRS10 data. We included in the comparison analysis the widely applied ACPC web-tool. The assessment took into account 27 clear-sky Landsat 8 scenes over a radiosonde station in Southern Brazil.

The obtained results showed that for the three atmospheric parameters (atmospheric transmittance, upwelling radiance, and downwelling radiance) the ACPC provided parameters in better agreement with those calculated using the radiosondes. The second-lowest uncertainties occurred when using CFSv2 profiles. No significant statistical differences were found between the parameters from the two WRF grids. None of the profile sources outperformed the others in all case days analyzed. The overall metrics of the WRF profiles were influenced by some cases with large errors. In general, the assessed profile sources better represented the atmosphere in dry and cold conditions. With respect to the retrieved LST values, using those calculated with the SBPA profile as reference, 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. Both WRF grids and CFSv2 generated a positive LST bias, while ACPC generated a negative bias. On balance, all the profile sources presented relatively good results in estimating the LST.

From the above findings, we conclude the following.

• 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.
• The results reinforce the validity and feasibility of ACPC, which is free of charge.
• Even though the overall statistical metrics for WRF profiles were inferior, their results were satisfactory for the estimation of both atmospheric parameters and LST values.

Despite the fact that some studies used the WRF model to simulate the skin temperature [118,119], to the best of our knowledge, this paper is the first attempt to apply the WRF model to aid the atmospheric correction of thermal remote sensing data. Its use showed potential, and our findings encourage further validations. Our study adds to the background of studies combing TIR satellite images and high-resolution NWP models.

Finally, in spite of the discoveries in this study, some limitations to be considered in future work are worth mentioning. There were no in situ LST data available to perform validation of the retrieved LST. The uncertainties in the SBPA radiosonde profiles were not analyzed; they were directly assumed as reference data. The study results refer to one area only, and therefore additional research is currently ongoing to extrapolate the conclusions to other environments.