Assessment of High-Resolution LST Derived From the Synergy of Sentinel-2 and Sentinel-3 in Agricultural Areas

This work explores the potential of obtaining high-resolution thermal infrared (TIR) data provided by the Sentinel-2 (S2) & Sentinel-3 (S3) constellation in a typical semiarid agricultural environment. Maps of land surface temperature (LST) with 10–20 m spatial resolution were obtained from the synergy S2–S3 in the Barrax test site in Spain, for a set of 14 different dates in the summers of 2018–2019. Ground measurements of LST transects covering a variety of croplands and surface conditions were used for a ground validation of the disaggregation approaches. A cross validation of the LST products was also conducted using Landsat-8/TIRS images. Two recent approaches exploiting the linkages between shortwave and thermal data were adapted and tested, with differences in the inputs, the physical-mathematical framework, or the treatment of the LST residuals, and two options for the original 1 km S3 LST data were considered. Despite the large range of temperatures registered (295–330 K), differences with observed values resulted in an average RMSE < 3.0 K and a negligible systematic deviation, showing good results even in small fields ∼1 ha. Results confirm the need for appropriate adjustment techniques of the LST residuals obtained to better capture the low temperature conditions. The systematic overestimations introduced by the use of the operational sea and land surface temperature radiometer L2 LST product, and the limitations associated with certain irrigation management are discussed. Results in this work offer a solution to the lack of high-resolution satellite TIR data, and provide new opportunities for LST applications in agricultural areas.

Landsat program has released an operational LST product as part of its Collection 2 [16], [17].Rescaled 30 m LST data are now available globally from the TIR sensors on board Landsat 4 to 9, but the 16-day revisit time (8-day revisit time with two Landsat satellites in tandem) might not be sufficient for certain goals.On the other hand, satellites such as Sentinel-3 or Suomi-NPP provide daily thermal observations but with much coarser spatial resolution.The Copernicus high spatio-temporal resolution land surface temperature monitoring mission [18] is expected to fulfill the spatio-temporal requirements of many agricultural applications, with a pixel size <50 m and a revisit time of less than 3 days.Other planned TIR missions with high spatio-temporal resolution are the thermal infrared imaging satellite for high-resolution natural resource assessment [19] and the surface biology and geology [20].But in the meantime, downscaling methods are contributing to bridge this gap by providing LST at a fine resolution [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32].
The combination S2-S3 accomplishes the spatio-temporal resolution necessary for agricultural applications mentioned above.The visible and near infrared (VNIR) information provided by the multispectral instrument (MSI) on board S2 can be used to characterize TIR subpixel heterogeneity within the S3 1 km pixel.Several papers have been recently published on this topic [39], [40], [41].The assessment of these techniques has been conducted by cross validation using LST products derived from Landsat or ASTER observations.However, robust validations based on ground-level LST measurements are scarce, due to the difficulty in compiling a ground LST dataset covering a variety of surface conditions within the same scene.Validation of satellite derived LST is essential to quantify their precision and accuracy.
This article follows-up previous research by Bisquert et al. [21] and Sánchez et al. [32], testing downscaling techniques in the semi-arid agricultural area of Barrax test site (Albacete, Spain).The combination of large parcels with small land holdings, irrigated and rainfed plots, and a wide variety of crops, and thus coverages, confers this study site unique features for LST validation purposes [17], [32], [42], [43], [44].
Two different disaggregation approaches were evaluated in this work.The disaggregation approach reformulated by [32] for the combination MODIS/S2, was now adapted to the tandem S3-S2, using the 1 km LST L2 product from SLSTR as input, but maintaining the same algorithm framework to derive 10 m LST maps (hereafter referred as S&G_10 m).The other algorithm to be explored was proposed by [39], based on a machine learning formulation to generate 20 m LST maps (hereafter referred as G&N_20 m).
Operational S3 LST product has been validated in a variety of sites [33], [36], [45], [46], [47] and an LST overestimation has been observed by these authors, attributed to wrongly assigned biomes linked to wrong coefficients in the SW algorithm.With the aim of exploring the effect of these uncertainties in the original S3 LST product, a self-derived SW algorithm was introduced and applied, as part of this work.
Ground LST data from two experimental campaigns carried out in the summers of 2018 and 2019 were used for the model assessment.A set of 14 Sentinel-3 scenes, with concurrent ground transects and favorable sky conditions, was selected for this work.A detailed analysis of the performance of both models was conducted, with particular attention to the following.
1) Overall performance assessment based on quantitative model validation using the ground LST measurements as a basis.2) Impact of the field size on the assessment results.
3) Model constraints based on surface condition changes due to irrigation.4) Effect of the local uncertainties in the original S3 LST L2 product used as input.The rest of this article is organized as follows.Section II describes the study site, and the ground and satellite data.A brief description of the two disaggregation approaches is also included in this section.Section III shows the results of the ground-based validation and the cross validation with L8-TIRS scenes.The analysis and discussion of the results is included in Section IV.Finally, Section V concludes this article.

A. Study Site and Measurement
This work was carried out in a semi-arid agricultural area of Southeastern Spain.The selection of the study site was based on the following research needs: variety of crops and land covers within a confined area of several km 2 , presence of fields of different extension (from <1 to >20 ha) as well as different irrigation systems (sprinkler and drip) together with rainfed plots, easy walking access into the cropfields for the transects collection, and prevalence of cloud-free weather conditions during spring-summer to maximize the availability of satellite scenes.The Barrax test site (39°03 N, 2°06 W, 720 m a.s.l.) accomplishes the above criteria, and this work focuses on ground data collected in "Las Tiesas" experimental farm and surroundings (see Fig. 1).
Suitability of Barrax area as a cal/val test site is well-known, as it has been used in many international campaigns [42], [43], [44], [48], and it is particularly attractive to evaluate the performance of disaggregation techniques [21], [32].
This research focuses on the data gathered during the summers of 2018-2019.Ground measurements of radiometric land surface temperature (LST g ) were taken concurrent to Sentinel-3 overpasses, covering several experimental plots per date.Up to 15 different locations/plots were selected, representative of a wide variety of crops and surface conditions.A set of four hand-held infrared radiometers (IRTs) Apogee MI-210 were available for this experiment, in addition to a multispectral radiometer CIMEL Electronique CE 312-2 [49].All radiometers were calibrated before and after the experimental campaigns.Ground temperatures were measured with an average accuracy of ±0.2 and ±0.1 K using Apogees and CIMEL, respectively.Transects were conducted carrying back and forth the IRTs, nadir pointing the surface from a height of 1.5-2.0m, at a rate of 5-10 registers/min, and covering as large area as possible within 10 min centered on the S3 overpass time.This resulted in coverage of a 60 × 60 m 2 area, representative of a grid of 3 × 3 (20 m) or 6 × 6 (10 m) Sentinel-2 pixels.All measurements were corrected for atmospheric and emissivity effects [17], [32], [50].Downwelling sky radiances were measured as part of each transect, and a land surface emissivity dataset gathered by the authors by combining the temperature and emissivity separation procedure [51] and CIMEL data were used to tailor correction parameters to each site/date.This was essential since surface/cover conditions in a specific plot/site may change drastically depending on the phenological phase of a crop.For instance, emissivity spectrum of a cereal cropfield differs a lot from the green phase to senescence or once tilled.LST g data were obtained as the average value of the 50-100 punctual IRT measurements registered in every plot/site within the 10 min frame centered in the S3 overpass.The standard deviation of these measurements was used to represent the spatio-temporal variability of the ground data.Regarding the woody crops, such as the vineyards or almond orchards, temperatures for soil and canopy components covering both shadow and sunlit portions were registered and carefully weighted to derive a target temperature [17].Note the majority of the selected plots are provided with a water supply system (sprinkler or drip), and are then subject to irrigation events during the experiment, except plots labeled as "2.2.Bare soil" and "3.1.Rainfed barley."This allows us to assess the impact of irrigation on the performance of the S3-S2 disaggregation.This will bring new insights in the feasibility of disaggregated LSTs for agricultural applications and help to understand its limitations.Furthermore, four of the sites have an extension <4 ha, with 2 of them limited to 1 ha.These are challenging conditions in which potential of S3-S2 disaggregation will be also explored.

B. Satellite Images
Sentinel-2A & 2B and Sentinel-3A & 3B images were used as inputs in this work.MSI on board Sentinel-2 has 13 spectral bands (ranging 443-2202 nm) with a spatial resolution of 10, 20, or 60 m.The combination of both S2A & S2B satellites guarantee a geometric revisit time of 5 days globally.L2A product containing bottom-of-atmosphere (BOA) reflectance was used in this work.These L2A products are generated from the associated L1C products, and are available over Europe since March 2018.These L2A BOA products are composed of 100x100 km 2 tiles in cartographic geometry.
Thermal data are acquired by the SLSTR on board Sentinel-3 satellites.SLSTR is provided with two TIR bands, with a spatial resolution of 1 km and a revisit time of 1-2 days.The Copernicus SLSTR Level-2 LST product was used as a first approach.In order to minimize the viewing angle (VA) effect, only scenes with VA<25°were selected.This product contains LST on the wide 1 km grid, as well as associated normalized difference vegetation index (NDVI) values for each gridded pixel, in one annotation dataset.Reflectance data from the shortwave OLCI instrument, also onboard S3, were used in this work.Full resolution Level-2 surface reflectances are provided at approximately 300 m pixel size.In particular, the Level-2 VGK NDVI product was used.
Sentinel-2 L2A and Sentinel-3 L2A products were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu,last accessed 16/10/2023).Single pairs of concurrent S2 and S3 images, or within ±1-day timing difference, were selected to run the S&G_10 m approach (see dates in Table I) whereas all the cloud free S2 scenes falling within 10 days before and after the S3 overpass were selected for processing in the G&N_20 m algorithm.MSI bands 4 and 8 were used to compose an NDVI image at a spatial resolution of 10 m in the framework of the S&G_10 m approach, whereas bands 1 to 7, 8 A, 11 and 12 were involved in the G&N_20 m algorithm.Fig. 2 shows an example of the spatial distribution of the NDVI over the study site.Note the wide range in NDVI values available in the area during the experiment.
Three Landsat-8/TIRS images, concurrent with S3 overpasses, were also used to perform a cross validation of the LST results.Images were downloaded from the USGS archives (https://earthexplorer.usgs.gov).Note the original spatial resolution of the L8 TIR bands is 100 m, although LST L2 product is provided in Collection 2 with a 30 m pixel size.These L8 scenes were used as a reference for the distributed assessment and analysis of the disaggregation techniques.

C. Downscaling Approaches
Both of the downscaling approaches rely on physical relationships which exist between shortwave and thermal observations.The strongest of those is the correlation between the amount of vegetation (observable using shortwave reflectance or indices such as NDVI) and the LST.In general, the denser the vegetation the cooler the surface, due to evaporative cooling (latent heat flux) but also due to increase in surface roughness (which increases turbulence and thus enhances heat transfer from the surface), shadowing, and other secondary effects.By exploiting reflectance in specific shortwave bands additional relations could be derived.For example, S2 bands 11 and 12 are around the shortwave infrared water absorption bands and therefore can provide information on leaf and top-soil water content which inversely impacts the LST (the higher water content, the higher evapotranspiration and the lower the LST).Similarly, S2 bands 5, 6, and 7 cover the red-edge parts of the spectrum which can provide information on the chlorophyll and thus on the potential photosynthetic activity of the plants, which again is inversely proportional to LST (the higher the photosynthetic activity the higher transpiration and the lower the LST).Those relationships are very variable in both space and time and therefore are specifically derived for each pair of shortwave and thermal observations.Sections below provide more details of the two downscaling approaches.
1) S&G_10 m: This approach is based on the sharpening method explored by Agam et al. [24], that was adapted by Bisquert et al. [11] to be applied to the combination MODIS/Landsat and further tailored to be applied to the tandem MODIS/S2 to derive 10 m LST maps [32].This algorithm has been now adapted to the combination S3-S2.
Briefly, NDVI values for the 10 m S2, NDVI_S2_10 m (calculated from reflectance values in the Red (R4, 665 nm), and NIR (R8, 842 nm) bands), and the 300 m S3, NDVI_S3, are aggregated to equivalent 1000 m S3 TIR pixel.Differences between S2 and S3 VNIR data due to spectral resolution, atmospheric correction, viewing angle or pixel footprint were corrected through a normalization extracted from the 1000 m NDVI, then applied to 10 m S2 NDVI (NDVI_S3_10 m).The 1000 m coarse spatial resolution required a previous selection of "pure" pixels for the NDVI-LST adjustment.This selection was based on a confidence value calculated from the comparison between NDVI_S3 and aggregated NDVI_S3_10 m.This confidence value was computed as the ratio between the standard deviation from the 3 × 3 pixels NDVI_S3 belonging to each 1000 m pixel, and its mean value, as suggested by Kustas et al. [2].Pixels with confidence values within the lowest quartile were selected in this step.A linear regression was established between NDVI_S3 and LST_S3 at 1000 m, using data from those "pure" pixels (1), and then applied to the NDVI_S2_10 m values to obtain a prime estimate of 10 m LST [LST_PRIME_10 m, (2)].An updated residual (Res) correction was proposed by Sánchez et al. [32] to account for the local conditions, and to correct the possible deviations produced by the NDVI-LST equation.This residue is calculated as the difference between the original and predicted LST at a coarse resolution (3), and further smoothed based on a linearization between the residual Res and the NDVI_S3 itself from 1000 m data (4).This linear relationship between the residue and the NDVI was then applied to 10 m NDVI_S2_10 m (5), see Fig. 3. Finally, 10 m LST values were obtained by adding this residual (Res) to original 10 m LST_PRIME_10 m data (6).This updated protocol to derive the residue values was shown effective in [32], reducing the LST deviation, mainly in small size fields surrounded by a different cover Res _1000m = LST_S3 _1000m − LST_PRIME _1000m (3) LST _10m = LST_PRIME _10m +Res _10m (6) where a 0 , a 1 , b 0 , and b 1 are the adjusted parameters.The data aggregation is a key step in the disaggregation procedures.The aggregation of the VNIR bands and NDVI was carried out by averaging the values of all high-resolution pixels within an equivalent low-resolution pixel.
Following paper [26], the aggregation of the TIR band was done through the Stefan-Boltzmann law with the assumption of similar emissivity values for adjacent pixels 2) G&N_20 m: This method is based on paper [26] and was previously applied by Guzinski et al. [39], [52] to sharpen thermal data to be used as input into evapotranspiration models.Each S3 scene is matched with a mosaic of all the S2 scenes acquired at most ten days before and after the S3 acquisition, and the regression model used for sharpening is derived specifically This method, implemented as an open-source Python application (https://github.com/radosuav/pyDMS,last accessed: October 16, 2023), relies on a machine-learning algorithm that derives a statistical relationship between high-resolution variables aggregated to low-resolution pixels and the low-resolution variables that need to be sharpened (i.e., LST).This derived relationship is then applied to the high-resolution data.An ensemble of decision-tree regressors is used in the current implementation of the algorithm, summarized as follows.
Briefly, the atmospherically corrected Sentinel-2 optical data with a spatial resolution of 20 m is resampled to match the pixel sampling of the SLSTR sensor (1 km).Concurrently, the SRTM DEM is used to derive slope and aspect maps which, together with S3 overpass time, are used to estimate the sun beam incident angle of a flat tilted surface.The DEM and the solar incidence angle maps are also resampled to the SLSTR resolution.Note this specific step will have a minor impact on this work due to the flat terrain of the Barrax area.A multivariate regression model is then trained with the three resampled datasets used as predictors and the LST used as the dependent variable.The selection of training samples is performed automatically by selecting the most homogeneous samples at 1 km scale by means of the coefficient of variation (CV) of all the high-resolution pixels falling within one low-resolution pixel, selecting 80% of pixels with the lowest CV where subscript i represents the spectral band, n is the total number of spectral bands, ρ and σ are the mean and the standard deviation of the fine resolution reflectances within the coarse resolution pixel, respectively.The regression model is based on a bagging ensemble [53] of decision regression trees.The decision trees are additionally modified such that all samples within a regression tree leaf node are fitted with a multivariate linear model, as proposed by Gao et al. [26].The regression models are trained on the whole S2 tile (100 km by 100 km) as well as on tiles of 30x30 S3 pixels.Once trained, they are also applied on the whole scene and on each window.The bias between the predicted high-resolution LST pixels aggregated to the low-resolution and the original low-resolution SLSTR LST is calculated and the outputs of the whole-scene and local-window regressions are combined based on a weight inversely proportional to the bias [26]: where p i represents the prediction from the local or global model, w ic is the weight at coarse resolution, and the subscript c represents coarse resolution.The prediction (local or global) with the smaller residual at a given coarse resolution pixel represents a better prediction and thus is weighted higher using the reversed squared residual as follows: where r ic is the residual for the local and global model prediction at coarse resolution.Finally, the LST predicted by the regression model is corrected by comparing the emitted longwave radiance of the sharpened fine LST versus the original coarse LST.A bias-corrected LST is therefore recalculated by adding an offset all fine scale pixels falling within coarse scale pixel in order to remove any residual bias.This is done to ensure the conservation of energy between the two thermal images with different spatial resolutions [26].The output of the sharpening is a 20 m representation of the LST.

D. Split-Window Algorithm to Reduce Uncertainties in Coarse Resolution S3 LST Inputs
The SW algorithm used in this work is an adaptation to small viewing angles of the algorithm introduced by Pérez-Planells et al. [36] for Sentinel-3.This is based on the algorithm for the spinning enhanced visible and infrared sensor on board METEOSAT Second Generation 2 presented in [54], and can be summarized as follows: where T stands for LST and BT 8 and BT 9 are at-sensor brightness temperatures in K for the SLSTR channels at 11 and 12 μm, respectively; ε is the mean emissivity for those SLSTR channels and the difference between them is noted as Δε; θ is the sensor viewing angle; and W is defined as the total water vapor column provided by LST L2 product.Emissivities were obtained using the NDVI-based approach proposed by Valor and Caselles [55].
The coefficients in (11) were obtained from regression analyses between LST -BT 8 and BT 8 -BT 9 , and using the blackbody approach (ε = 1 and Δε = 0) following [56] procedure and the CLAR database.This radiosoundings database is composed of 382 clear-sky atmospheric profiles distributed over all latitudes, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and covering a wide range of vapor water content (0-7 cm) and temperatures (253-313 K).In this work, Gaussian angles from 0°to 40°were chosen to generate the dataset for training the SW algorithm.A total of 10 696 different cases were used to obtain the coefficients in (11).

A. Ground Validation
Maps of disaggregated LST from both models, S&G_10 m and G&N_20 m, were obtained for all available dates (see Fig. 5).Results of the comparison between disaggregated LST and ground measurements (LST g ) are shown Fig. 6 and Table II.Values for the disaggregated LST from both models, S&G_10 m and G&N_20 m, correspond to a 3 × 3 pixel average centered in each measurement transect, whereas the central pixel was extracted for the S3_L2 LST data.Error bars correspond to the standard deviation of the 3 × 3 pixel average, in the y-axis, and to the standard deviation of the ground measurements, in the x-axis, representing the spatial variability of the transects and the temporal homogeneity of the 10 min frame concurrent to the S3 overpass time for each scene.Temperatures range between 295 and 330 K, with the lowest LSTs corresponding to high NDVI conditions and the largest LSTs to bare soil conditions.Following paper [32], data after any sprinkler irrigation event registered were excluded from this first analysis to avoid disruption in the expected thermal behavior.
Results for >70 sample data show a slight overall overestimation of 1.4 and 1.8 K for S&G_10 m and G&N_20 m, respectively, with an RMSE of 3.3 and 4.0 K.The improvement in comparison with original S3 LST is clear, with an RMSE of 7.2 K and an overestimation of 3.2 K affecting land surface temperatures extracted from the 1x1 km 2 pixels.Note estimation errors below 1.5 K were reported by Pérez-Planells et al. [36] or Yang et al. [45] using the operational S3 LST product in favorable and homogeneous surface conditions.
These statistics agree with the linear regression parameters, with a correlation coefficient of 0.87 and 0.83 for S&G_10 m and G&N_20 m, respectively.Fig. 6 shows model overestimations mainly occurring for lower LST conditions.Furthermore, results were grouped in terms of surface conditions using a crop classification as a basis (see Fig. 7).Differences between disaggregated and measured LST remain between ±4.0 K for most classes.An overall lower scatter within each crop type was observed using G&N_20 m, in comparison to S&G_10 m.This might be a positive effect of considering all the S2 bands in the G&N_20 m approach versus the single NDVI information in the S&G_10 m algorithm.However, the overestimation with G&N_20 m was stressed for potato, fescue and corn croplands, in comparison with S&G_10 m.The reason might be linked to the capacity to discriminate LSTs values in small parcels where the difference between the 10 and 20 m pixel size may have an effect.

B. Distributed Assessment
The ground-based vicarious calibration was further completed with a cross calibration of the disaggregated scenes using Landsat 8 Collection 2 LST product as a basis for the assessment.Three golden dates, July 17, 2018, July 24, 2018, and June 18, 2019 were selected with concurrent L8 and S3 overpasses, with no more than 15-20 min acquisition delay.Disaggregated S3-S2 LST maps were resampled to the spatial resolution of the L8 LST data.Pixel-to-pixel differences were calculated at 30 m spatial resolution for a selected subset of 5 × 5 km 2 covering our study site (see Fig. 8).Note the significant improvement in the parcels discrimination at the disaggregated 10-20 m resolution compared to the original L8 LST pixel size.
Differences range between ±10 K, with a Gaussian/normal distribution centered within ±1.3 K and an RMSE below 3.0 K for the three dates and both, S&G_10 m and G&N_20 m approaches.These results are in agreement with the vicarious calibration above.And again, no significant differences are observed between the performance of both models.
However, discrepancies between disaggregated S3-S2 LSTs and L8 product data are not equally distributed throughout the scene.Differences are particularly pronounced for those areas covered by sprinkler-irrigated crops (summer irrigated crops), Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.for which disaggregation approaches still fail at capturing the cooling down effect subsequent to a recent irrigation event.This will be further discussed as follows.
Despite some limitations, findings from these point-based and distributed assessments are promising for the application of both, the prime NDVI-LST linkage, and the more sophisticated decision tree regressions involving all reflectance bands, to the synergy S2-S3 for the monitoring of high-resolution LST in croplands.

A. Impact of Extreme Temperatures
Fig. 5 illustrates the failure of Sentinel 3 in capturing the full thermal variability within the scene, due to spatial aggregation within each 1 km pixel.This may have an impact on the capacity of the disaggregation approaches to resolve extreme thermal conditions.This is stressed for highly vegetated fields, typically characterized by low LST values, and not so evident for the opposite bare soil, and high LST values, conditions.Results from original S3 plotted in Fig. 7 reinforce this idea, showing the largest overestimations for the festuca, potato and corn fields, and the main underestimation for the bare soil plot.This dichotomy is also maintained for the G&N_20 m outputs, although not so evident for the S&G_10 m approach.Fig. 9 plots these deviations as a function of the absolute values of LST g registered for both approaches.The higher uncertainties occur for the cropfields mentioned above where the temperatures are always below 310 K.
Note the extension of these critical parcels is < 4 ha.However, the field size itself is not the main reason for these LST deviations since the 10-20 m spatial resolution of S2 is able to resolve details in fields < 0.5 ha.The vineyard is a good example since its intermediate temperatures, ranging from 305 to 320 K, are reproduced with an average estimation error lower than 3.0 K, in agreement with the overall error for the full dataset, despite its field size < 4 ha.
For a more in-depth analysis of thermal discrimination capability, the statistics for the model performance were recalculated  III show the results from both models improve after excluding those three parcels, with an overestimation and RMSE of 0.9 and 3.0 K, respectively, for S&G_10 m, and 0.6 and 2.6 K for G&N_20 m.Note a better agreement with ground measurements is achieved by G&N_20 m algorithm compared to S&G_10 m under these conditions, with the original S3_L2 LST providing an average estimation error of 4.2 K.These findings show the potential of including all VNIR 20 m S2 bands as inputs in the disaggregation algorithm (G&N_20 m) in comparison to the sole use of NDVI (S&G_10 m).However, focusing on those data from corn, potato, and fescue, results for G&N_20 m show a bias of 5.7 K and RMSE of 6.6 K, whereas a bias of 2.9 K and RMSE of 3.8 K are obtained with S&G_10 m.Differences with original S3_L2 products are also very significant with a bias of 11.7 K and RMSE of 12.7 K.These findings indicate S&G_10 m performs better for highly vegetated surface conditions, which might be a consequence of the additional residual correction step, proposed by Sánchez et al. [32] and implemented in this work to correct the deviations produced by the NDVI-LST equation.In any case, improvement versus the original TIR spatial resolution of S3 is clear, and derived LST are robust and feasible under a variety of surface conditions.Moreover, these results indicate the pixel size (10-20 m) is not determinant for the accuracy of the disaggregation approach, even in parcels <1 ha.

B. Constraints Due to Crop Irrigation
Table III includes results for those data affected by a sprinkler irrigation event within 12 h prior the S3 overpass.The overestimation reaches ∼10 K for disaggregated LST with no significant differences between both algorithms, nor with the original S3_L2 values.The larger time scale for the S2 images used to match each S3 scene in G&N_20 m (± 10 days) versus S&G_10 m (± 1 day) algorithms does not result in any significant degradation at this point.
A land cover classification map, derived from S2 multispectral imagery, was used to identify and extract these parcels (summer irrigated crops), shown in Fig. 10.Overestimation increases with the temporal proximity of the irrigation event, and peaks for those areas with concurrent water supply during the L8 overpass.Note a good performance of the disaggregated products is also possible in these areas, when the plots were not recently irrigated, or the temporary cooling effects disappear.This is the case in most fields for 06/18/2019 (see Fig. 10).
Nevertheless, disaggregation approaches can capture LST fairly well over drip-irrigated cropfields.As an example, Fig. 10 includes LST-difference results filtered for almond orchard areas, using the classification map.This comprises both dripirrigated and rainfed almond orchards, with differences ranging ±3.0 K, with no appreciable differences between both treatments.These results are promising for rainfed, but also drip-irrigated, woody crops such as almonds or pistachio orchards, since S3-S2 disaggregated LST can help in increasing the frequency of daily ET estimates through surface energy balance modeling in these crops.

C. Effects of the Uncertainties in the S3 LST L2 Product
In a recent work, [36] showed a systematic uncertainty of 1.5 K and a precision of 1.2 K in the S3 LST L2 product, using ground measurements over a rice paddy site, with long tradition as a cal/val site [50], [54], [57], [58], [59].Similar values were reported by other authors.For instance, Li et al. [47] found an average overestimation of 1.7 K, together with a RMSE of 3.4 K, using the official S3 LST product against in situ LSTs from six stations in the HiWATER experiment in Northwest China.In [34], an uncertainty of 1.9 K was obtained at a dessert area in Namibia, with a bias of 1.8 K.An average error of 1.4 K was reported by Zheng et al. [46], with a bias of 1.1 K. Yang et al. [45] obtained an RMSE of 2.4 K, with a bias of 1.6 K in the same dessert of Namibia.This LST overestimation is attributed to wrongly assigned biomes linked to wrong coefficients in the SW algorithm.For instance, Pérez-Planells et al. [36] showed accuracy improvement in LST using local-adjusted coefficients, reducing the systematic deviation to 0.4 K and the RSD to 1.1 K. Other works obtained a bias of -0.8 K with a correctly assigned biome, as rainfed croplands in this case [34].Yang et al. [45] Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.reported an RMSE of 0.7 K, with a bias of 0.4 K over a large water body, correctly classified as water body.Recently, Li et al. [47] analyzed the performance of various LST retrieval algorithms and highlighted the retrieval errors due to the effect of the land cover type misclassification on the surface emissivity values.According to these findings, the operational L2 LST product for S3 is affected by atmospheric and emissivity uncertainties in the implemented algorithm, and adjustment to local conditions may improve LST accuracy.For this reason, a self-derived split-window algorithm (introduced in Section II-D) was applied in this work, based on coefficients calculated from the CLAR database [56].S&G_10 m disaggregation approach was run using these SW LSTs as a basis, to evaluate the effect in LST accuracy.Overall results improve, but particularly for small parcels (<4 ha), as shown in Fig. 11 and Tables II-III.Systematic deviation results negligible now, and RMSE decreases to ±2.9 K. Focusing on smallest plots <1 ha, RMSE remains now stable with a slight underestimation of 0.7 K.
These results highlight the needs for accurate S3 LST products as inputs in the disaggregation algorithms to avoid propagation of undesirable uncertainties into final high resolution LST estimates.

V. CONCLUSION
Data fusion techniques between S3 and S2 platforms have an evident potential to bridge the existing gap in high-resolution and high-frequency TIR imagery.Whether disaggregated LST data from this tandem can serve as input to models providing ET information for agricultural applications is not fully answered in the literature, due to the lack of robust LST ground datasets in these ecosystems.The heterogeneity and variability in surface conditions of the Barrax experimental site provides a unique opportunity to address this point.Cross validation using Landsat or ASTER TIR images is useful for a distributed assessment, but feasibility in agricultural applications needed to be tested with a comprehensive dataset of ground LST transects under a variety of field crops, irrigation systems, and field size patterns.This assessment concludes high-resolution LST can be extracted from the S2-S3 tandem in croplands with an overall accuracy of ±2.6 K and negligible systematic deviation.Nevertheless, strengths and limitations of this study need to be exposed as follows.
1) Based on the results of this work, including all VNIR 20 m S2 bands as inputs in the disaggregation scheme can benefit the LST estimates compared to the integration of the sole NDVI information.However, the better performance of the S&G_10 m approach for high vegetated surface conditions indicates the importance of implementing a residual correction step in the process to account for deviations produced by the reflectances-LST relationship.2) This study indicates the pixel size (10-20 m) is not determinant for the accuracy of the disaggregation approach.Results of LST disaggregation for ∼1 ha plots, with an RMSE below 3.0 K, show the robustness of the S&G_10 m approach also under these challenging conditions.3) S3/S2 disaggregation fails in the areas affected by a recent sprinkler irrigation event prior the S3 overpass, due to the cooling down of the surface.The natural relation between LST and surface reflectance response cannot capture the decrease in surface temperature registered by ground thermal radiometers, that may reach >10-15 K.The overestimation effect disappears when the irrigation is applied 1-2 days before.4) Results are positive for rainfed, but also drip-irrigated cropfields where the cooling down effect is minor, and the disaggregation approaches perform well.This is the case of woody crops such as vineyards, almonds, or pistachio orchards.These findings are of particular interest in semiarid areas with scarcity of water resources, and a growing preference for these crops, since S3-S2 disaggregated LST will help in monitoring their water status and water needs through surface energy balance modeling, with no critical restrictions due to the parcel size.5) Accuracy in S2-S3 disaggregated LST certainly improves by tuning the coarse resolution S3 LST data.An adjustment of the coefficients in the S3 LST algorithm to the local conditions is highly recommended to minimize effects of wrong atmospheric or emissivity characterization in final high resolution LST estimates.In this work, as much as 1 K overestimation in S3 LST in our area is removed by using a self-derived SW algorithm.Although further research could contribute by replacing reflectances or NDVI by other indices accounting for the leaf water-or chlorophyl content, or even including some ancillary soil moisture data, findings in this work are encouraging for the use of these disaggregation methodologies in agroecosystems.Results are promising since the S2 10-20 m pixel size, together with the daily revisit frequency of S3 satellites, can fulfill the LST input requirements in a variety of hydrological, climatological, and agricultural applications.However, this work shows evidence that the agroecosystem needs will not be fully satisfied until high spatio-temporal TIR missions are orbiting and operationally capturing the biophysical variability.

Fig. 1 .
Fig. 1.Overview of the study site.Ground measurement points are located over a S2 false color composition corresponding to date May 16, 2019.Labels for the different study fields are explained in the adjacent Table.White circles illustrate the 1 km diameter areas enclosing the study sites.

Fig. 2 .
Fig. 2. Examples of S2 NDVI images to illustrate the phenology evolution in the study site during the 2018/19 experiment.Dates are stamped in the upper right corner of each scene.

Fig. 3 .
Fig. 3. Flowchart of the S&G_10 m downscaling methodology, including the different processing steps, inputs and outputs.Variable descriptions are included in the text.

Fig. 4 .
Fig. 4. Flowchart of the G&N_20 m downscaling methodology, including the different processing steps, inputs and outputs.Variable descriptions are included in the text.

Fig. 5 .
Fig. 5. Disaggregated LST from S&G_10 m (right column) and G&N_20 m (central column), together with the original S3_SLSTR (left column) products.Examples correspond to a subset of 5 × 5 km 2 covering our study site, and dates: July 17, 2018 (top row), July 24, 2018 (central row), June 18, 2019 (bottom row).Labels for the corresponding LST products and dates are stamped in the lower left and upper right corners, respectively, of each scene.

Fig. 6 .
Fig. 6.Linear regression between disaggregated LST (both, S&G_10 m and G&N_20 m) and ground-measured values (LST g ).Dashed line represents the 1:1 agreement.Error bars correspond to the standard deviation of the 3 × 3 pixel average, in the y axis, and to the standard deviation of the ground measurements, in the x-axis, representing the spatial variability of the transects and the temporal homogeneity of the 10 min frame concurrent to the S3 overpass time for each scene.

Fig. 6
Fig.6plots the linear regressions for the three LST products.Values for the disaggregated LST from both models, S&G_10 m and G&N_20 m, correspond to a 3 × 3 pixel average centered in each measurement transect, whereas the central pixel was extracted for the S3_L2 LST data.Error bars correspond to the standard deviation of the 3 × 3 pixel average, in the y-axis, and to the standard deviation of the ground measurements, in the x-axis, representing the spatial variability of the transects and the temporal homogeneity of the 10 min frame concurrent to the S3 overpass time for each scene.Temperatures range between 295 and 330 K, with the lowest LSTs corresponding to high NDVI conditions and the largest LSTs to bare soil conditions.Following paper[32], data after any sprinkler irrigation event registered were excluded from this first analysis to avoid disruption in the expected thermal behavior.Results for >70 sample data show a slight overall overestimation of 1.4 and 1.8 K for S&G_10 m and G&N_20 m,

Fig. 7 .
Fig. 7. Box plot of the differences between disaggregation results and ground LST measurements for each land use: S&G_10 m (top) and G&N_2 0m (central).Differences between original S3-L2 LST and ground measurements are also included (bottom).Outliers are labeled with a nonfilled circle and X represents the mean value.Shaded areas correspond to plots with extension <4 ha.

Fig. 8 .
Fig. 8. Differences between disaggregated LST and L8 LST product: G&N_20 m (central column) and S&G_10 m (right column) for dates: July 17, 2018 (upper row), July 24, 2018 (central row), June 18, 2019 (bottom row).Landsat 8 LST products are also shown (left row).Labels for the corresponding LST products and dates are stamped in the lower left and upper right corners, respectively, of each scene.

Fig. 9 .
Fig. 9. Differences between modelled and measured LST versus the absolute values of LST g .Filled marks correspond to the reduced dataset (corn, potato, and fescue).Shaded areas indicate the average estimation errors for the full dataset.

Fig. 10 .
Fig. 10.Examples of the differences between disaggregated LST and L8 LST product, G&N_20 m (left column) and S&G_10 m (right column), filtered for sprinkler-irrigated areas (top row images) and almond orchards (bottom row images).Labels for the corresponding LST products and dates are stamped in the lower left and upper right corners, respectively, of each scene.

TABLE I LIST
OF SENTINEL-3 SCENES USED IN THIS STUDY AND SUMMARY OF THE METEOROLOGICAL CONDITIONS IN THE AREA AT THE OVERPASS TIME

TABLE III STATISTICS
OF THE DIFFERENCES BETWEEN DISAGGREGATED LST OR S3_LST, AND GROUND-MEASURED LST DATA, BY CONSTRAINING EXTREME CONDITIONS by isolating those data (corn, potato and fescue) from the full dataset.Statistics listed in Table