NOAA-AVHRR Orbital Drift Correction: Validating Methods Using MSG-SEVIRI Data as a Benchmark Dataset

: National Oceanic and Atmospheric Administration–Advanced Very High Resolution Radiometer (NOAA-AVHRR) data provides the possibility to build the longest Land Surface Temperature (LST) dataset to date, starting in 1981 up to the present. However, due to the orbital drift of the NOAA platforms, no LST dataset is available before 2000 and the arrival of newer platforms. Although numerous methods have been developed to correct this orbital drift effect on the LST, a lack of validation has prevented their application. This is the gap we bridge here by using the 15 min temporal resolution of Meteosat Second Generation–Spinning Enhanced Visible and Infra-Red Imager (MSG-SEVIRI) data to simulate drifted and reference LST time series. We then use these time series to validate an orbital drift correction method based on solar zenith angle (SZA) anomalies that we presented in a previous work (C1), as well as two variations of this approach (C0 and C2). Our results show that the C0 method performs better than the two others, although its overall bias absolute value ranges up to 1 K, while standard deviation values remain around 3 K. This is veriﬁed for most land covers, for all NOAA platforms, and these statistics remain mostly stable with noise on SZA time series (from 0 ◦ to ± 10 ◦ ). With this study, we show that orbital drift correction methods can be thoroughly validated and that such validation should aim toward bias absolute values below 0.1 K and standard deviation values around 1.4 K at coarse spatial resolution. To validate other orbital drift correction approaches, the drifted and reference time series used in this work are freely available for download from the ﬁrst author’s webpage. This will be the ﬁrst step toward the building of an orbital-drift-corrected long-term LST dataset.


Introduction
Land surface temperature is a key parameter for the monitoring, modeling, and forecasting of water and carbon cycles in our context of climate change. Temporally coherent measurements of this parameter are available since the year 2000, thanks, for example, to the moderate-resolution imaging spectroradiometer (MODIS) onboard the Aqua and Terra platforms [1,2]. Before 2000, the Advanced Very High Resolution Radiometer (AVHRR), onboard the successive National Oceanic and Atmospheric Administration (NOAA) platforms, is the only instrument (see Table 1) to allow for the estimation of this parameter [3], although platform characteristics prevent their use because of temporal incoherence [4]. Therefore, it is-at the moment-impossible to further relate observed worldwide air temperature increases to observed sea and land surface temperatures, as already carried out recently from the MODIS archive [5,6].
Actually, none of the NOAA platforms was designed to remain on its nominal orbit, so the platforms' orbit started to drift immediately after their positioning [4]. This resulted in a progressively later acquisition of data, introducing a spurious decreasing trend in afternoon land surface temperature (LST) time series throughout the lifetime of each NOAA satellite. This is the so-called orbital drift effect, which has been identified in [4]. It has little influence for surfaces with high heat capacity, especially water, although for arid or cultivated areas [7], this drift can introduce in the data a decreasing trend of various degrees, masking easily the underlying trend of expected temperature increase due to climate change. Several methods have been presented to correct this orbital drift effect [8][9][10][11][12][13][14][15]. These can be classified into three groups: sun-target-sensor geometry modeling, LST daily cycle reconstruction, and statistical methods. The first group uses vegetation structure models to estimate the proportion of shaded and lit areas to estimate each pixel LST for each NOAA overpass [11]. This approach has been used to build a daily record of NOAA-14 AVHRR land surface temperature over Africa [16]. The second group relies on the reconstruction of the LST daily cycle, which allows, once identified, to retrieve the LST that would have been measured at a given hour [8,10,15]. The reconstruction of the LST daily cycle is the key aspect of this approach and can be carried out from modeling [10], interpolation from several measurements through the same day [8], or the local neighborhood [15], under the assumption of local similarity of surface characteristics. The third group relies on time series analysis, using the drifted LST time series to isolate the drift from a proxy signal, usually the solar zenith angle (SZA), to represent the image acquisition time [9,[12][13][14].
Despite these numerous correction methods, no global orbital-drift-corrected LST dataset has yet been compiled, although some datasets allow for LST estimation, such as the Long-Term Data Record (LTDR [17]), as shown in [7]. The main reason for this resides in the difficulty in validating this corrected dataset. Indeed, validating such a dataset would be immensely difficult due to the absence of in situ LST data for a long period covering at least part of the 1980s and the 1990s, over different locations of our planet, with different land covers. However, validating the orbital drift correction may be feasible, even if only a few attempts exist in the literature, relying mainly on time series statistics improvement [10,14].
As a matter of fact, recent geostationary platforms allow for an estimation of the LST at short time intervals (less than 30 min) and at a spatial resolution close to the one of most global datasets (a few kilometers). Actually, we can select the geostationaryderived LST acquisitions at orbital drift simulated times to build a drifted LST time series, apply an orbital drift correction to this time series, and compare this corrected time series to a constant-time geostationary-derived LST time series (for example at 13:30 UTC -Coordinated Universal Time). This is the approach chosen here through the use of LST data estimated from Meteosat Second Generation (MSG) Spinning Enhanced Visible and Infrared Imager (SEVIRI) data [18,19].
More specifically, the aim of this paper is to present an approach for the validation of orbital drift correction methods, applied to methods based on the one proposed by the authors in a previous work [14]. More specifically, this work will answer the following questions:

•
How should we build a reference and validation LST dataset for orbital drift effect correction assessment (Sections 2 and 3)? • How do different methods for orbital drift effect correction fare with regard to time series stability requirements (Section 4)?

Data
Since 2007, the Global Change Unit of the University of Valencia is receiving MSG-SEVIRI data continuously, through a direct broadcast high-resolution picture transmission (HRPT) system, which consists of a parabolic dish, a personal computer with the hardware and software to decode L-band data, and a set of storage devices to save all the received SEVIRI data. Received data amount to 96 images per day for each channel due to the 15 min temporal resolution of SEVIRI, resulting in storage needs of 1TB per month. To build our validation and reference time series corresponding to the different NOAA platforms, we used data from 2013 to 2019. Missing data, due mainly to maintenance operations, were retrieved from the EUMETSAT (European Organisation for the Exploitation of Meteorological Satellites) data center as high-rate information transmission (HRIT) raw files and then processed following the same scheme as the data retrieved from our station.
MSG-SEVIRI provides data every 15 min, in 12 spectral channels: 4 in the visible and near-infrared wavelengths (0.4-1.6 µm) and 8 in the infrared wavelengths (3.9-13.4 µm). These bands are termed high-resolution visible (HRV), VIS06, VIS08, IR016, IR039, WV062, WV073, IR087, IR097, IR108, IR120, and IR134. Bands WV6.2 and WV7.3 are located in a water vapor (WV) absorption region, whereas band IR9.7 is located in the ozone absorption region and band IR134 in a CO 2 absorption region. The other bands are located in regions with high atmospheric transmissivity (atmospheric windows), and in particular, bands IR108 and IR120 are located in the split-window region, and are used for surface and cloud observation. Spatial resolution of the SEVIRI instrument is 1 km at the nadir for the highresolution visible (HRV) channel and 3 km for the other channels. Over southern Europe, the typical extension of a pixel is 4-by-5 km. Due to the geostationary character of the MSG platform, only half the planet is imaged, centered on the geographical coordinates (0 • ,0 • ). The whole processing scheme of the data received by our MSG-SEVIRI station is described in [18], while the LST and SST (Sea Surface Temperature) algorithms are presented in [19]. In summary, this process consists of the inversion of Planck's equation to retrieve brightness temperatures for IR108 and IR120 bands, from which surface temperatures are estimated through the split-window approach, through a previous estimation of emissivities and total atmospheric water vapor [19].
We also downloaded MSG-SEVIRI cloud masks from the EUMETSAT archive for the years 2013-2019. These masks are built following the methodology developed by [20]. The masks are provided for pixels with a viewing zenithal angle below 80 • , whose values 0, 1, 2, and 3 correspond, respectively, to clear over sea, clear over land, cloudy, and no data. In this work, we only used LST estimations corresponding to cloud mask values of 1. Cloud mask values are unavailable for pixels with a viewing zenith angle above 73 • before 2016.
For visualization purposes, and to save computer memory as well as processing time, we selected all BEnchmark Land Multisite ANalysis and Intercomparison of Products (BEL-MANIP [21]) sites included in the MSG-SEVIRI disk, provided their viewing zenithal angle was below 73 • (177 in total). The locations of these pixels can be viewed in Figure 1, along with their corresponding International Geosphere Biosphere Programme (IGBP) class.

Methods
We present here the methodology used to carry out our objectives, which are split in three parts: building of the validation time series, presentation of the orbital drift effect correction methods, and statistical methods used for method assessment. First, we used the equations provided in [22], which allow one to estimate NOAA platforms' equatorial crossing times for their activity period. Here, we used the parameters provided at https://www.star.nesdis.noaa.gov/socd/sst/3s/ (accessed on 1 March 2021) for the following platforms: NOAA-7, NOAA-9, NOAA-11, NOAA-14, NOAA-16, NOAA-18, and NOAA-19 (termed hereafter N07, N09, N11, N14, N16, N18, and N19). For any coordinates (latitude, longitude) and date, these equations give the hour at which a given NOAA platform has crossed the equator (local time). These equations were retrieved from orbital parameters retrieved from NOAA ephemeris and Two-Line Elements (TLE) files. We provide here the equatorial crossing time (EXT) equation (Equation (1)): where t corresponds to time (in days) since the platform launch; η 0 is the average EXT (in hours); α 1 , ω 1 , and ϕ 1 are the amplitude (in hours), frequency (in days −1 ), and phase (in radians) of the first harmonic, respectively; and α 2 , ω 2 , and ϕ 2 are the amplitude, frequency, and phase of the second harmonic of the fit, respectively (for details, see [22]). For the up-to-date values of the fit coefficients (η 0 , α 1 , ω 1 , ϕ 1 , α 2 , ω 2 , ϕ 2 ), see https: //www.star.nesdis.noaa.gov/socd/sst/3s/ (accessed on 1 March 2021). Figure 2 presents the equator crossing times for the N07, N09, N11, N14, N16, N18, and N19 platforms. From this equatorial crossing hour (solar time), we need to estimate the UTC hour of NOAA platform overpass for our selected pixel (T(t,x,y)), which is given by Equation (2): where Φ(x,y) and λ(x,y) are, respectively, the latitude and longitude (in degrees) of our given pixel (x,y) and δ is the inclination of the satellite orbit (99 • for NOAA platforms). Since our MSG-SEVIRI LST dataset is limited to the years 2007-2019, and to respect possible seasonal effects on the time series, we chose to match the last available date of each NOAA platform (as it appears in the LTDRV4 dataset: https://ltdr.modaps.eosdis.nasa. gov/cgi-bin/ltdr/ltdrPage.cgi, accessed on 1 March 2021) with the corresponding day of the year (DOY) during 2019. As a consequence, all LST measurements of our validation series (drifted and reference) lie within the period 2013-2019.
Then, for each BELMANIP site, we identified the two MSG-SEVIRI images immediately preceding and following this UTC hour of pixel overpass and their corresponding LST values and, by linear interpolation, retrieved the MSG-SEVIRI LST for the overpass time. If any of these two observations is labeled as cloudy in the corresponding EUMETSAT MSG-CMA product, we then labeled the retrieved LST value as cloudy. By repeating this for all dates of each platform activity period, we obtained, for each BELMANIP pixel, the drifted time series. For the reference time series, we proceeded similarly, starting from the reference local time of 13:30 (solar time) and then proceeding as above from Equation (2) onward. Thus, we obtained, for each BELMANIP pixel lying within the observational disk of MSG-SEVIRI, the reference and drifted LST time series for the 7 NOAA platforms mentioned above. These time series are freely available for download from the first author's webpage (www.uv.es/juy/resources.htm, accessed on 1 March 2021).
The methods we selected for validation are based on a previous work of ours, described in [14]. This approach is pixel based and uses anomalies of different parameters (SZA, LST) to identify the part of the LST time series that is due to the orbital drift effect. These anomalies are calculated as the difference of the parameter time series and its average year, estimated as the sequence of LST averages over each DOY: for 1 January, the average of all LST values retrieved on 1 January over the lifetime of the chosen NOAA platform; for 2 January, the average of all LST values retrieved on 2 January; and so on until 31 December. For leap years, LST values retrieved on 29 February were added to the ones of 28 February to calculate the average LST for the latter DOY.
The first approach we tested is the one described in [14]. This approach relies, for each NOAA platform, on a quadratic fit of the SZA anomaly against time. Then, the LST anomaly is fitted against this SZA fit and a linear function of time. Finally, the part of the LST anomaly explained by the SZA fit is subtracted from the LST time series. For a thorough description of this correction method, please refer to [14]. The second approach we tested is a variation of the previous one, where the LST anomaly is fitted against the SZA fit alone (no dependence on time). The third approach is also a variation of the first method, this time with the LST anomaly fitted against both the SZA fit and a quadratic function of time. In the following sections, the LST time series retrieved after these corrections are, respectively, labeled LST C1 , LST C0 , and LST C2 , with the numerical index corresponding to the degree of the polynomial fit of the LST anomaly against time, while the drifted and reference time series are, respectively, labeled LSTD and LSTR. As for the corresponding correction methods, we termed them C1, C0, and C2, respectively.
To validate the orbital drift effect correction methods, we estimated for each satellite and each BELMANIP site displayed in Figure 1 the bias, standard deviation (STDV), and root-mean-square error (RMSE) between the corrected time series (LST C0 , LST C1 , or LST C2 ) and the reference time series (LSTR). To check whether the correction does improve the time series, we also estimated these statistics between the drifted (LSTD) and the reference (LSTR) time series. We then averaged these statistics for each IGBP land cover. We processed similarly to validate the orbital drift corrections for each NOAA platform. To assess the influence of the MSG-SEVIRI viewing angle, we also estimated these statistics Finally, since the AVHRR instrument has an across-track viewing angle of ±55 • , which introduces variation in sun-target-sensor geometry, we repeated the process described above with the addition of white noise of different amplitudes (0 • to 10 • by steps of 1 • ) to the drifted SZA time series. This white noise is created by a different uniformly random process for each BELMANIP pixel and NOAA platform, and its influence on the 3 orbital drift correction methods are assessed with the above-mentioned statistics. Figure 3 presents the drifted and reference time series for the N07 platform for a selected desert pixel located in eastern Libya, labeled BELMANIP site 208. The top graph of Figure 3 shows the LST annual behavior of a desert pixel, with LST values ranging between 290 and 330 K, with a slightly decreasing trend for the drifted time series (continuous) not observed for the reference time series (dotted). The difference between both time series reaches up to 8 K at the end of the simulated activity period (4 years). This is shown in the Figure 3 bottom graph, which presents the difference between drifted, C0, C1, and C2 corrected time series, on the one hand, and the reference time series for the same site and the same platform, respectively, in red, green, blue and black, on the other hand. We observed that all correction methods improve on the drifted time series, with C1 and C2 correction methods performing similarly and C0 slightly better. Note that the initial difference is not null, due to the difference between the N07 initial equatorial crossing time (around 14:30 solar time) and our chosen reference (13:30 solar time). For this specific site and platform, statistics for the C0, C1, C2, and drifted time series are presented in Table 2.

Results
These statistics show a clear improvement after correction, the C0 approach performing better than the two other ones, which show similar results. However, in the best case, the remaining bias is still -1.3 K, around the accuracy of the LST estimation, while the initial bias was −4.1 K, leading to differences between drifted and reference time series up to 12 K for the N07 simulated activity period (Figure 3, top graph). Therefore, lower corrected bias values are needed for adequate correction of the orbital drift for this site.   Table 3 presents the overall statistics for all sites and all platforms. Once again, C0 correction performs better with regard to the RMSE and standard deviation statistics, although C1 correction performs better with regard to bias removal. As in the case of site 208 for the N07 platform, the standard deviation decreases slightly after correction, while bias values are closer to zero. Since these three characteristics are complementary, we focused hereafter on the bias and standard deviation statistics, the aim of a perfect orbital drift correction being to reach null bias and to decrease slightly the standard deviation. Of course, reaching a null standard deviation is illusory due to short-term variations in the land surface temperature (changes in wind power and direction, alternating presence of sun and clouds), which cannot be eliminated by orbital drift correction methods. Table 3. Overall statistics for all selected BELMANIP sites and all NOAA platforms between drift-ed and reference time series as well as for each orbital drift correction method (C0, C1, and C2). STDV stands for standard deviation.  Figure 4 presents the bias statistics for the three orbital drift correction methods tested here (C0, C1, and C2) as well as the bias of the drifted time series (LST D ) for each BELMANIP site. This figure shows that the highest (absolute values above 2 K) original bias values (from the drifted time series) are located in all land covers, except the forested areas of Amazonia, Central Africa, and northern Europe. These biases generally decrease (in absolute value) for all corrections, although C0 correction seems to perform better than the two other correction approaches. More specifically, bias values for C0 correction remain highly negative in southern Africa, the Sahel, and western Brazil. These highly negative values are also found for C1 correction, although a few highly positive biases are also evidenced for several sites. These highly positive bias values are not found in the case of C2 correction, even though they remain higher (in absolute values) than for C0 correction.  C2 (b, c, and d, respectively). Figure 5 presents the standard deviation values for all sites for the three orbital drift correction methods (C0, C1, and C2) as well as the bias of the drifted time series (LST D ). With regard to the standard deviation between drifted and reference time series, the lowest values (below 3 K) are found in forested areas, while the highest values are found in arid and semi-arid areas. Note that in most cases, the standard deviations remain below 5 K, except for a few sites, probably with high spatial heterogeneity. All corrections decrease slightly the standard deviation values, except for C1 correction in a few sites, where the standard deviation remains the same or increases slightly. Figure 6 summarizes these results as a histogram of the bias and standard deviation values obtained for the drifted time series as well as the three corrected time series. In Figure 6, and all the following graphs, the graph for the drifted time series (LST D ) is plotted as black diamonds, while the graph for the corrected time series (LST C0 , LST C1 , and LST C2 ) are, respectively, plotted as red triangles, green squares, and blue crosses. While standard deviation values after correction are only slightly lower than for the drifted time series, there is a clear difference with regard to bias statistics. C1 correction presents the histogram best centered around 0, although C0 and C2 corrections have a narrower distribution, C0 correction bias being closer to 0 than C2 correction.  C2 (b, c, and d, respectively). When we look at these statistics by IGBP land cover (Figure 7), we see that the orbital drift effect is highly land cover dependent, with higher effects for bias (below −2 K) and standard deviation (above 3 K) for given classes corresponding to shrublands, savannas, grasslands, and deserts (namely classes 5, 6, 7, 8, 9, 10, and 16), as observed previously in the literature [7]. Figure 7 confirms that C0 correction performs better than the others, with the resulting bias above −1 K and standard deviations slightly lower than the drifted ones. As for C1 correction, bias values are better centered around 0, as commented above, while standard deviation values are slightly higher than for the drifted time series for many classes (1, 2, 4, 5, 6, 7, 13, and 14). This, combined with positive bias values, shows an over-correction of the orbital drift. As for C2 correction, its bias and standard deviation values are located halfway between the drifted and C0 correction values.  We observe that the drifted time series present higher absolute values (for both statistics) for the first NOAA platforms, namely N07, N09, and N11, while the lowest absolute values are found for N18 and N19 platforms. As before, C0 correction performs better than C1 and C2 corrections with regard to standard deviation, while C0 and C1 corrections perform similarly with regard to bias, except for N07 and N11. C2 correction performs worse than C1 correction for most platforms with regard to bias, although it performs slightly better with regard to standard deviation. Once again, C1 correction tends to over-correct the orbital drift, especially for the N11 platform.  Figure 9 shows how the introduction of noise in the drifted SZA time series influences the correction methods with regard to bias and standard deviation. As a matter of fact, C0 and C1 corrections are highly stable for both statistics, while C1 correction tends to over-correct the orbital drift effect more with increased noise (increase in both bias and standard deviation).  Figure 10 presents the percentage of cloud-free data averaged over all platforms for each BELMANIP site. As expected, the availability of data is higher (up to 75%) for arid areas and lower (down to 15%) for temperate and equatorial areas. Table 4 presents the correlation values between the drifted and corrected statistics retrieved above and the cloud-free percentage. These correlations remain low (absolute values below 0.24 for corrected time series), pointing toward a low influence of cloud contamination on orbital drift correction performance.

Discussion
We have seen that although C0 correction performs better than the other two corrections and improves over the drifted time series, the difference between the corrected time series and the reference time series remains high, with bias absolute values up to 1 K, depending on the considered platform or land cover. While C1 correction sometimes present better bias values than C0 correction, it sometimes over-corrects the orbital drift effect, resulting in increased standard deviation. This happens when the fit of the LST anomaly against time and the SZA anomaly are compensated one by the other, as mentioned in [14]. However, the LST data used here come from the SEVIRI sensor, while the AVHRR LST data, to which the correction will be applied, present a few differences with regard to LST characteristics.
The first difference is the spatial footprint: the SEVIRI LST pixel covers 3-by-3 km at the nadir, widening poleward, to 4-by-5 km, for example, over southern Europe. As for AVHRR, its footprint is of 1.1-by-1.1 km, although datasets usually provide the data at lower resolution, for example 0.05 • -by-0.05 • for LTDR [17], which translates roughly to a 5-by-5 km pixel. Since the order of magnitude of pixel extension is similar, and since land cover is usually heterogeneous at scales above 1 km [23], this should not hinder our validation scheme.
Another difference between SEVIRI and AVHRR data is the absence of changes in sensor-viewing geometry. While the SEVIRI observation angle is fixed, the AVHRR scans the Earth below across its track, with the observation angle varying between ±55 • . This results first in varying AVHRR pixel dimensions and therefore varying pixel heterogeneity. This aspect may be partially accounted for during the reprojection step of dataset building. This varying observation angle results also in variations in the retrieved LST due to surface anisotropy, in addition to the varying amount of sunlit and shaded proportions of the observed surface [11]. Finally, the varying observation angle changes the atmospheric path of the signal, although some LST retrieval algorithms account for this aspect through a previous estimation of the total amount of atmospheric water vapor [3]. This aspect can be simulated through an added white noise on LST time series.
Moreover, using Equation (1) for each BELMANIP pixel assumes that this pixel is located at the nadir of each NOAA platform for each acquisition. This hypothesis is obviously not met, since the AVHRR has a swath of 2900 km. This introduces an incertitude of the order of ±51 min in the equator crossing time estimate, corresponding to one orbital period. This incertitude can be simulated by white noise on the equator crossing time estimates, or directly on SZA values, as carried out here. We saw (Figure 9) that such noise had no influence on the quality of the reconstruction for C0 and C2 methods, while the C1 method shows performs the worst with increased noise amplitude.
For this work, we used the activity periods for NOAA platforms as displayed in the LTDRV4 dataset. However, a newer version of this dataset has been recently made available (LTDRV5), with similar activity periods, for most satellites, except for N19, whose activity period is a lot longer, starting in 2009 up to today (almost 12 years). Figure 8 shows that almost no drift is evidenced for this platform (up to the end of 2015), although Figure 2 shows that N19 suffers a stronger orbital drift after 2015. Therefore, the statistics values for this platform will probably increase with a longer period activity in order to become comparable with the N11 platform. However, this platform, as well as N18, shows that when the orbital drift is controlled (by well-designed orbit choice and not-too-long activity period), the resulting bias between the drifted and reference time series is of the order of 0.1 K, while standard deviation values are around 1.4 K. These are therefore the bias and standard deviation values toward which correction methods should aim for all land covers and all platforms.
Finally, this validation is only a validation of the presented orbital drift effect correction methods. However, this does not guarantee a valid time series after reconstruction, because the successive AVHRRs were not placed on the same orbits, with nominal equator crossing times varying from 13:30 to 14:00 (solar time) for NOAA afternoon satellites. Therefore, a sensor intercalibration step is mandatory for the building of an orbital-drift-effect-corrected LST dataset. Additionally, comparison of the corrected LST time series with in situ data is recommended, although, as stated in the introduction, this is a difficult task, since in situ time series are scarce for earlier periods of the satellite era (1980s and 1990s). Moreover, the difference between the spatial extension of the in situ measurements and the AVHRR pixel size also makes the validation difficult. To illustrate this point, we can take a look at the only study (to our knowledge) that tried to validate the orbital drift effect correction method at the pixel level [15]. This validation resides in the comparison of the error between in situ LST data and the AVHRR LST before and after the orbital drift effect correction. This error drops from 3.5K to 3.0K after the correction (average RMSE estimated from their study's Figure 11 from [15]). While the improvement is real, the relative difference is low, even if, as mentioned by the authors, such high error values are common in large pixel LST validations [15]. However, our approach allows the estimation of the error generated by the orbital drift correction method alone, independently of other sources of error (detection noise, surface heterogeneity, etc.).
The validation and reference time series used here have been built by extracting from our MSG-SEVIRI dataset the LST data for BELMANIP pixels, along with the corresponding SZA and cloud mask values. These data are freely available for download from the first author's webpage (www.uv.es/juy/resources.html, accessed on 1 March 2021). These data are useful for the validation of any pixel-based orbital drift effect correction method and could be extended relatively easily (upon demand) to the correction methods based on the estimation of the LST daily cycle from local neighborhoods (such as [15]).

Conclusions
This study presents the first thorough validation for orbital drift effect reconstruction methods from LST time series. It has shown that even though the validation of the orbitaldrift-effect-corrected LST is difficult, the validation of the correction method can be carried out using MSG-SEVIRI geostationary data, with similar spatial resolution and 15 min temporal resolution. This validation has shown that the C0 method performs better than the other two correction methods, with bias absolute values up to 1 K, depending on the considered platform or land cover, and a standard deviation of 3 K. Additionally, we showed that an optimal orbital drift correction should aim toward 0.1 K bias and 1.4 K standard deviation values and that cloud contamination has a low influence on correction performance. By using the BELMANIP sites for this validation, we studied how these corrections perform for each land cover, with visible improvement for the land covers such as deserts, shrublands, savannas, and grasslands. Although the tested statistics show that there is room for improvement, the presented LST time series provide a benchmark for testing most of the existing orbital drift reconstruction methods. These time series can also be used to design correction methods by analyzing the relation between drifted and reference time series of both the SZA and the LST. However, such validation is only a first step toward building an orbital-drift-corrected LST dataset, which should be completed by in situ LST validation of the corrected time series for as many as possible representative sites. The reference and validation data used in this study are available freely from the first author's website (www.uv.es/juy/resources.html, accessed on 1 March 2021).