Downscaling Meteosat Land Surface Temperature over a Heterogeneous Landscape Using a Data Assimilation Approach

A wide range of environmental applications require the monitoring of land surface temperature (LST) at frequent intervals and fine spatial resolutions, but these conditions are not offered nowadays by the available space sensors. To overcome these shortcomings, LST downscaling methods have been developed to derive higher resolution LST from the available satellite data. This research concerns the application of a data assimilation (DA) downscaling approach, the genetic particle smoother (GPS), to disaggregate Meteosat 8 LST time series (3 km ˆ 5 km) at finer spatial resolutions. The methodology was applied over the Crau-Camargue region in Southeastern France for seven months in 2009. The evaluation of the downscaled LSTs has been performed at a moderate resolution using a set of coincident clear-sky MODIS LST images from Aqua and Terra platforms (1 km ˆ 1 km) and at a higher resolution using Landsat 7 data (60 m ˆ 60 m). The performance of the downscaling has been assessed in terms of reduction of the biases and the root mean square errors (RMSE) compared to prior model-simulated LSTs. The results showed that GPS allows downscaling the Meteosat LST product from 3 ˆ 5 km2 to 1 ˆ 1 km2 scales with a RMSE less than 2.7 K. Finer scale downscaling at Landsat 7 resolution showed larger errors (RMSE around 5 K) explained by land cover errors and inter-calibration issues between sensors. Further methodology improvements are finally suggested.


Introduction
Land surface temperature (LST) is a key variable for a wide range of environmental applications, including urban surface temperature mapping ( [1]), surface heat island monitoring ( [2]), and wildfire detection ( [3]), to cite but a few.More common applications involving LST are the monitoring of water and energy budgets, such as deriving surface energy budgets ( [4,5]), assessing drought severity ( [6]), vegetation physiology ( [7]), evapotranspiration ( [8,9]), and land surface model (LSM) calibration ( [10,11]).Nevertheless, reliable monitoring of such environmental applications requires LST measurements at frequent intervals and fine spatial resolutions to capture LST space/time variability.Despite the advances of remote sensing, measuring LSTs at both high spatial and temporal resolutions is not yet possible considering current space sensors.While high spatial resolution thermal data have limited temporal availability (e.g., a revisit period of 16 days for Landsat 7 is insufficient for agronomic applications considering that only clear sky observations are of use), coarser satellite thermal data are available at higher frequencies (e.g., 15 min temporal resolution for Meteosat 8).However, coarse resolution pixels are often subject to mixture effects that could be presented as a blending of multiple temperature components within the large thermal pixel ( [12]).This mixture effect is relatively more prevalent in spatially heterogeneous environments, such as agricultural ( [13,14]) and urban areas ( [1]).The main challenge is, therefore, to develop methods to assess LST at higher spatial and temporal resolutions.
A well-established method to derive higher resolution LSTs from other available data is disaggregation, also known as downscaling ( [12]).The basic tenet of downscaling is to produce finer resolution sub-pixel temperatures from the coarser resolution thermal data through sharpening or unmixing algorithms.Different approaches have been so far developed, ranging from statistical regression schemes ( [13,15,16]), to data assimilation (DA)-based ones [14,[17][18][19].Statistical regression-based approaches are built on the inverse relationship between LST and vegetation indices, such as the Normalized Difference Vegetation Index (NDVI) or the fractional vegetation cover (fcover), both increasing with the amount of biomass ( [13,16,20]).Nevertheless, LST variability cannot only be explained by vegetation amount and depends also on several other factors, such as photosynthetic activity of vegetation ( [21]), surface albedo, soil moisture, or emissivity conditions ( [22][23][24]).Furthermore, statistical regression-based techniques are hard to apply in an operational data production system because they rely on empirical relationships and stationarity assumptions, difficult to extend spatially.Moreover, different inputs may be important over different landscapes making difficult the choice of a prediction model which may vary from one location to another.Using a more robust LSM may give more flexibility for a global application of the downscaling scheme.DA-based techniques could then be a smart alternative to derive higher LST, as demonstrated in [14,17,19].The main idea of DA approaches is to constrain a physical modeling of the sub-pixel temperatures provided by the LSM with large scale observations based on a prior knowledge of vegetation class fractions inside the thermal pixel.In fact, under classical assumption of Gaussian error, unbiasedness, and linearity of observational operator, the optimal estimation is obtained as the minimizing state of a quadratic cost function.DA ensures the estimation of the "optimal" sub-pixel temperatures through the minimization of a cost function.For downscaling purposes, the cost function relies on the discrepancy between coarse resolution observations and the aggregated modeled sub-pixel LSTs (more details are available in [19]).In real applications, the non-gaussianity of the errors, which is due to the non-linearity of the dynamic and/or of the observational operator, leads to a sub-optimal estimation.An optimal solution can still be obtained from the non-linear filtering theory, where a practical implementation is known as particle filtering [25].For this DA algorithm, the optimal state is the average of an ensemble of states weighted according to their proximity to the observations following the Bayes rule.However, in presence of bias, the particle filter only provides a sub-optimal estimation.
In this work, a DA-based downscaling approach has been used to enhance the spatial resolution of Meteosat 8 SEVIRI LST time series over the Crau-Camargue region in Southeastern France.The DA approach is based on a genetic particle smoother (GPS) which has been implemented into the SETHYS LSM [10,11,26,27].Previous studies have also dealt with the problem of assimilating coarse Earth Observation data and the need to disaggregate these (e.g., [19]).A step further is presented here, by evaluating the proposed method over larger time/space scales in real conditions using remote sensing products at various spatial and temporal resolutions.

Methodology and Data
This study focuses on the application and evaluation of our GPS downscaling approach over the Crau-Camargue region in Southeastern France.The period under analysis spans from March 2009 to September 2009 and the original observations are LST time series provided by the EUMETSAT Satellite Application Facility for Land Surface Analysis (LSA-SAF), http://landsaf.ipma.pt/[28], derived from SEVIRI thermal infrared observations on-board Meteosat 8 (the first of the second-generation Meteosat satellites).The evaluation of the downscaling is performed quantitatively at both 1 km and 60 m spatial resolutions using Moderate Resolution Imaging Spectro-radiometer (MODIS) LST products and Landsat 7 thermal infrared data, respectively.The study region and the various datasets used are described further in Section 2.3.

GPS Downscaling Methodology
Our GPS downscaling method fully described in [19], is based on the assimilation of low-spatial/high-temporal resolution LST observations within a LSM able to simulate the temperature space/time variability at sub-pixel scales.It supposes that the spatial heterogeneity of LST is mainly related to land cover [15] and soil moisture variability, the latter depending mostly on precipitation and irrigation practices.The assimilation is then performed at the low spatial resolution (LSR) pixel scale, each pixel assumed independent.In other words, each low spatial resolution pixel is downscaled independently to its neighborhood without considering any spatial, nor temporal, correlation.Such assumption is legitimate since, at the low spatial resolution of the observations, the spatial and temporal variabilities are lower than at the sub-pixel scale due to the mixing effects.In this way, the downscaling problem turns to the estimation of the respective sub-pixel temperature corresponding to the main land cover classes present in the considered pixel, assuming the stationarity of the land cover fractions.

Data Assimilation Approach
The estimation of a state variable in DA requires the definition of two models: the dynamic model, which describes the dynamics of the state variables in space/time, and the observation model, which relates the observations to the state variables.In our case, the dynamic model is SETHYS ( [10]), which simulates the surface energy and water budgets (and, therefore, the surface fluxes, surface temperature, and soil moisture).The model is applied to the various land cover classes present in the LSR pixel to simulate prior values of the respective surface temperatures which are aggregated to assess the observable LSR temperature through the aggregation model corresponding to the observation model [19].The simulation of different sub-pixel temperatures in the SETHYS model is possible by varying some parameter values (the greatest influence on SETHYS surface temperature estimation) with respect to their physical ranges.Consequently, the assimilation of LSR observations allows a two-fold result: model parameter calibration and sub-pixel temperature estimation [19].
However, with LSR observations being a mixture of the different sub-pixel temperatures present within the pixel, different combinations of sub-pixel temperatures could be solutions of the downscaling problem.We are then faced with what is called an equifinality problem and ensemble DA methods present a good framework since the sub-pixel LSTs are estimated from an ensemble of candidate solutions ( [29,30]).
For our purpose, we work with the particle filter ( [31,32]), an ensemble method successfully applied to state and parameter estimation and downscaling issues ( [11,33]).A brief presentation of the methodology is given in Appendix A. The particle filter was implemented to calibrate a set of model parameters to minimize the misfit between the simulated temperatures and observations.In practice, the misfit is quantified through the computation of the cost function deduced from the log-likelihood.In our case, since the time constant of LST is short and the observations can be irregular (due to cloudiness), the misfit is calculated over a time window, or assimilation window, varying between one day to one week ( [11]).Therefore, a smoothing implementation of the particle filter has been proposed for that purpose.

Implementation of the GPS
In general, particle filtering includes two main steps: the analysis and the selection/re-sampling.The first step consists in weighting each particle (i.e., an ensemble member, realization of a state variable) according to the "closeness" of the simulations to the available observations (the closer, the higher the weight).For time-independent observational errors, the weight is first computed over the assimilation window as the product of the likelihood relative to each observation.A normalization step is performed so the final weights form a discrete probability density.In this case, we talk about smoothing rather than filtering, since the observations are considered over a time interval, which is the opposite of a sequential filter, where they are considered at a single time.The second step, which represents the "heart" of particle filtering, consists in rejecting low-weighted particles and replacing them with copies of the kept ones (re-sampling).In our case, a particle represents a set of SETHYS parameters relative to all the actual sub-pixels inside the LSR pixel (the most influential/sensitive parameters towards SETHYS surface temperature output) and the cost function is the "distance" between simulated LSR and observed ones.In other words, if we have a number of Q land cover classes (end-members) inside a LSR pixel and if, for each class "i", a number of K i parameters (among N param ) have been retained, the particle can be written as follows: x " pP 1  1 , P 2 1 , . . ., where x is the particle and P j i is one of the sensitive model parameters of class i.Once the particle is defined, the prior ensemble of particles can be generated after defining a probability distribution and a range of variation for each parameter.In our case, and in the absence of highly sampled in situ measurements, we opted for a non-informative distribution to explore the whole parameter space.A uniform distribution is a suitable choice since it allows exploring the particles (parameters) space with equal probability.The parameter intervals of variation were prescribed empirically according to previous work ([10,19,27]).Then, for each parameter set corresponding to the particle, SETHYS simulates prior estimations of the class temperatures, given the atmospheric forcing data and vegetation characterization.The aggregation step (the observation model), is performed finally to obtain an ensemble of simulated LSTs at the spatial resolution of the observations to be assimilated.The role of the particle smoother is to select the simulations which best fit the observations at LSR, the sub-pixel temperature solutions and parameter set values being the average over the posterior particle ensemble.The two-fold solution is then reached by the calibration of SETHYS parameters for each land cover class according to LSR available observations over the assimilation window and the downscaling of the LSR LST time series ensured with the estimation of the sub-pixel temperatures.We should note that the aggregation step implies the availability of an HSR land cover map, which provides the different land cover fractions inside each LSR thermal pixel.

Study Area and Forcing Data
The study area is located in Southeastern France (43.53 N 4.66 E, see Figure 1).The climate is typically Mediterranean, with strong winds and irregular rainfall (long dry periods in spring and summer).
This region is characterized by a wide range of land cover types and various irrigation practices.Micro-meteorological forcing data at half-hourly time steps, weekly vegetation measurements (Leaf Area Index (LAI) and vegetation height), and half-hourly soil moisture measurements necessary to initialize and force the SETHYS model, were provided from the "Crau-Camargue" database acquired by the EMMAH team over 2009 for several agricultural fields (the micro-meteorological stations specifically designed to sample the main land cover classes of the region, are presented with red location markers in Figure 1, upper panel).In this study, we focus on the spring and summer periods, from March 2009 to September 2009, because this period is crucial for crop temperature monitoring and water stress detection.Land cover was derived from a former image classification of Landsat 7 images (Courault, D., personal communication) in eight classes and validated against ground surveys.This classification was revised for this work in order to group the orchards and woods classes, which present the same soil and vegetation hydro-thermal parameters regarding the SETHYS model.After that step, seven classes were defined as follows: "Bare soil", "Water" (representing salt marches and swamplands), "Tree cover" (grouping orchards and woods), "Irrigated grasslands", "Summer crops (grouping sunflower, maize and sorghum), "Winter crops (grouping winter wheat and barley)", and "Rice" (see Figure 1, lower panel).
LAI and vegetation height, needed as inputs of SETHYS, were estimated for each of these classes, respectively, from NDVI data derived from the Landsat 7-TM reflectance data (at 30 m spatial resolution) and expert knowledge for some specific classes (tree cover, crops) over the studied period.Empirical relationships proposed by [34,35] and calibrated by [36] over part of the study area in 2006 have been used to infer LAI from Landsat 7 NDVI measurements over the year 2009.Daily values were finally interpolated after a careful merge with the available irregular field observations.Seasonal variations of the vegetation height were also estimated from experimental data and phased with the mean average seasonal variations of LAI per class.

Meteosat 8 SEVIRI Data
The Land Surface Analysis Satellite Application Facility (LSA-SAF) provided the low spatial resolution LST data that have been downscaled in this work.The LST product [27] was obtained with a generalized split-window algorithm ( [37]) from top of atmosphere brightness temperatures measured in the thermal infrared, namely in Meteosat 8 channels IR 10.8 μm and IR 12.0 μm.The LSA-SAF LST is produced at full Meteosat 8 spatial and temporal resolutions, with a 15 min Land cover was derived from a former image classification of Landsat 7 images (Courault, D., personal communication) in eight classes and validated against ground surveys.This classification was revised for this work in order to group the orchards and woods classes, which present the same soil and vegetation hydro-thermal parameters regarding the SETHYS model.After that step, seven classes were defined as follows: "Bare soil", "Water" (representing salt marches and swamplands), "Tree cover" (grouping orchards and woods), "Irrigated grasslands", "Summer crops (grouping sunflower, maize and sorghum), "Winter crops (grouping winter wheat and barley)", and "Rice" (see Figure 1, lower panel).
LAI and vegetation height, needed as inputs of SETHYS, were estimated for each of these classes, respectively, from NDVI data derived from the Landsat 7-TM reflectance data (at 30 m spatial resolution) and expert knowledge for some specific classes (tree cover, crops) over the studied period.Empirical relationships proposed by [34,35] and calibrated by [36] over part of the study area in 2006 have been used to infer LAI from Landsat 7 NDVI measurements over the year 2009.Daily values were finally interpolated after a careful merge with the available irregular field observations.Seasonal variations of the vegetation height were also estimated from experimental data and phased with the mean average seasonal variations of LAI per class.

Meteosat 8 SEVIRI Data
The Land Surface Analysis Satellite Application Facility (LSA-SAF) provided the low spatial resolution LST data that have been downscaled in this work.The LST product [27] was obtained with a generalized split-window algorithm ( [37]) from top of atmosphere brightness temperatures measured in the thermal infrared, namely in Meteosat 8 channels IR 10.8 µm and IR 12.0 µm.The LSA-SAF LST is produced at full Meteosat 8 spatial and temporal resolutions, with a 15 min sampling interval and a spatial resolution of about 3 km ˆ5 km on our region (3 km ˆ3 km at the sub-satellite point).The product is provided over land within the Meteosat 8 disk in clear sky conditions along with uncertainty estimation.An automatic quality control (QC) is performed on LST data, and the quality information is provided on a pixel-by-pixel basis.The LST confidence level was defined based on the following parameters: viewing angle, atmospheric characteristics (mostly column water vapor), and land surface emissivity (EM).The three considered levels of confidence (above nominal, nominal and below nominal) correspond to estimated uncertainties of LST values (respectively, less than 1 K, between 1 and 2 K, and above 2 K) and with an uncertainty threshold of 4 K.More details are provided in the product user manual available at http://landsaf.ipma.pt/algorithms.This uncertainty dataset was used to generate the observation error statistics required by the GPS.The LST time series processed spans from March to October 2009.

MODIS LST Data
The MODIS sensor onboard the Terra and Aqua platforms, offers a view of the Earth's surface within 36 spectral bands (0.4 µm, 14.4 µm).The MODIS (Terra/Aqua) Land Surface Temperature and Emissivity (LST/E) level 3 product (MOD11A1/MYD11A1) is a daily tile-based and gridded (sinusoidal projection) product at 1 km spatial resolution.The equatorial passing times for Terra and Aqua are, respectively, 10:30 a.m./p.m. and 1:30 a.m./p.m. leading to a minimum number of four observations per day [38] at the equator (more images can be acquired in the northern hemisphere because of the polar orbit).In this work MYD11A1 and MOD11A1 Collection 5 LST/E products were acquired over the study area for the same observation period (from 1 March 2009 to 30 September 2009) to evaluate the downscaling results at 1 km spatial resolution.More details are provided at http://modis-land.gsfc.nasa.gov/temp.html.For simplicity, we will use in the following term 'MODIS' to refer to 'MYD/MOD11A1 Collection 5 1 LST/E products.
As the MODIS LST product can merge different orbits, the acquisition time and viewing angles can be different from one pixel to the next.Thus, a pre-processing step of co-registration with Meteosat has been applied in order to compare observations acquired at similar time and to mask remaining cloudy pixels.For that purpose, only MODIS LST data with acquisition time-lag less than 10 min towards Meteosat data were kept, leading to 23 clear-sky MODIS LST dates kept for the evaluation process.

Landsat 7 Data
The Enhanced Thematic Mapper Plus (ETM+) on board Landsat 7 satellite acquires images of the Earth since July 1999 with a revisit frequency of 16 days.Landsat ETM+ has eight spectral bands (six bands in the visible, one thermal band, and one panchromatic band).An ETM+ scene has an instantaneous field of view (IFOV) of 30 m ˆ30 m in bands 1-5 and 7, while the thermal band 6 has an IFOV of 60 m ˆ60 m.The equatorial crossing time is nominally 10:00 a.m.˘15 min.In collaboration with the French National Center for Space studies (CNES), clear sky images have been extracted over 2009 from the USGS website (http://edcsns17.cr.usgs.gov/NewEarthExplorer/)over the study area and georeferenced according to the Lambert 3 projection.Atmospheric corrections have been performed to estimate LST from top of the atmosphere radiances using the MODTRAN5 radiative transfer model ( [39]) and comparison with ground measurements showed good agreement (see [40]).More details on the Landsat 7 data processing are available at http://www.cesbio.ups-tlse.fr/multitemp.The acquisition time over the study area is 10:20 A.M. and the seven following dates have been made available: 21  As with MODIS, we will use the term "Landsat 7" to refer to the ETM+ data used in the following.

Downscaling Application and Results
The GPS was implemented on the study region and applied to the Meteosat time series at the pixel scale to optimize within the pixel class temperatures.We first present the GPS particle definition chosen and how the ensemble simulations were performed.The downscaling experiment settings are then presented.

Particle Definition
As already noted, the particle definition (i.e., the calibration parameters) is very crucial since the downscaling solution relies on LST sensitivity towards the model parameters.In an early step of this work, the most sensitive SETHYS parameters were identified from a set of 22 parameters using a variance-based sensitivity analysis (Sobol sensitivity analysis, [41]) performed for each land cover type of the study area.The parameter intervals of the sensitivity analysis have been chosen based on the same intervals prescribed empirically according to previous works ( [10,19,27]).The results exhibited a maximum of five sensitive parameters for the considered land cover classes which are listed in Table 1.These parameters were used to define the particles in the downscaling procedure.It is important to point out that the prior temperature of the water class (mostly saltmarshes and swamps) has been assessed by considering this class as wet, bare soil forced to saturated soil moisture at each time step, therefore evaporating at a potential rate.

Downscaling Experiments
Table 2 describes the experimental settings of the different downscaling scenarios performed.Various sensitivity experiments to the observation and model errors, as well as the number and timing of the observations, were performed [19].The time period considered during the day (morning, afternoon, night) was also tested and proved to be significant because of the dynamic characteristics of the model and observation errors [19].For example, model errors appear to be larger at night due to the difficulties to correctly simulate the turbulence in more stable conditions.Therefore, we propose to limit the assimilation window to daytime (between 7:45 a.m. and 4:00 p.m., local hour) as already discussed in previous works [11,19].The impact of the atmospheric forcing uncertainty was also assessed by perturbing input data.We present here the results obtained for four experiments.Experiments 1 and 2 consider a daily assimilation window length and explore the impact of a perturbation on the input forcing data.In Experiment 1, the air temperature (which is the most influential atmospheric variable on LST modeling in SETHYS [11]) was perturbed using a first-order auto-regressive filter with an error variance equal to 1 K, i.e., this noise was additive and time-correlated to the air temperature.It increases when the air temperature trend is positive, and decreases when it is negative, and its value depends on the amplitude of the air temperature.Experiments 3 and 4 examine, respectively, the impact of the lengthening of the assimilation window to three and seven days, respectively.These time periods have been chosen following previous works [10,11], showing the importance of considering several consecutive days to calibrate the soil hydraulic parameters controlling the soil evaporation processes.It is also important to note that, for all experiments when more than 90% of particles are rejected at the end of the analysis step, all particles are re-sampled using a uniform distribution over the parameter space to avoid filter collapse problems.For each of these experiments, the performances of the GPS downscaling approach have been finally assessed at 1 km scale against the MODIS LST product for the 23 available dates and at a 60 m scale against Landsat 7 LST for the seven available dates over 2009.For these comparisons, pseudo-images have been reconstructed at the respective scales (1 km and 60 m) by aggregation of the weighted class temperatures using the land cover map at 30 m resolution.They were compared to the prior estimations provided by SETHYS.These prior solutions were calculated based on a "reference set" of model parameters based on expert knowledge of the average values of each considered vegetation and surface parameters.

GPS Downscaling LST Evaluation at 1 km Scale
Figure 2 presents the downscaling results for DOY 240 (28th of August at 11:00 UTC over a highly heterogeneous segment (22 km ˆ22 km) of the study area (red box plotted over the land cover map in Figure 1), when a MODIS image was available.The assimilated SEVIRI observations and the downscaled LST are compared to the MODIS product at 1 km resolution.The figure shows that the main spatial structures of LST are well captured by the GPS filter in coherence with the land cover input data.The largest temperatures are observed for the bare soil surfaces and the lowest for the tree cover areas.The MODIS image appears fuzzier because of co-registration errors and pixel overlapping, unlike the GPS reconstructed image based on the discrete land cover map.More quantitatively, Table 3 presents the statistics calculated between the GPS estimations and the observations, on the whole set of MODIS data (23 images) and for the four downscaling experiments.The prior model values, as well as the efficiencies, are also indicated to assess the contribution of the GPS.The efficiency (Eff) is defined as 1 minus the ratio of the posterior to prior RMSEs and is expressed as a percent.The maximum value is, therefore, equal to 1 and positive values indicate the benefits of the DA to lessen prior error statistics.First, the results show that the prior model simulations present a negative bias of −2.3 K (model cooler than the MODIS observations), that has been partly corrected with the data assimilation.The root-mean-square errors (RMSE) are also much improved in all of the experiments by more than 1 K.For the best case (daily assimilation window and noised input forcing), a value of 2.7 K is obtained, showing a significant improvement compared to the prior estimations at 4 K.A slight improvement of 0.1 K is obtained compared to the other experiments where no uncertainty on the forcing data has been prescribed.The efficiencies are all positive and range between 27.5% and 32.5%).This result is quite important because it proves the robustness of the downscaling method towards forcing uncertainties.It also shows that such perturbations may have a positive impact on the downscaling results when the model variability is poor.In fact, the later point can be explained by the fact that perturbing the input parameters enlarges the spread of the ensemble simulations and improves the finding of the sub-optimal solutions by increasing the degrees of freedom and consequently More quantitatively, Table 3 presents the statistics calculated between the GPS estimations and the observations, on the whole set of MODIS data (23 images) and for the four downscaling experiments.The prior model values, as well as the efficiencies, are also indicated to assess the contribution of the GPS.The efficiency (E ff ) is defined as 1 minus the ratio of the posterior to prior RMSEs and is expressed as a percent.The maximum value is, therefore, equal to 1 and positive values indicate the benefits of the DA to lessen prior error statistics.First, the results show that the prior model simulations present a negative bias of ´2.3 K (model cooler than the MODIS observations), that has been partly corrected with the data assimilation.The root-mean-square errors (RMSE) are also much improved in all of the experiments by more than 1 K.For the best case (daily assimilation window and noised input forcing), a value of 2.7 K is obtained, showing a significant improvement compared to the prior estimations at 4 K.A slight improvement of 0.1 K is obtained compared to the other experiments where no uncertainty on the forcing data has been prescribed.The efficiencies are all positive and range between 27.5% and 32.5%).This result is quite important because it proves the robustness of the downscaling method towards forcing uncertainties.It also shows that such perturbations may have a positive impact on the downscaling results when the model variability is poor.In fact, the later point can be explained by the fact that perturbing the input parameters enlarges the spread of the ensemble simulations and improves the finding of the sub-optimal solutions by increasing the degrees of freedom and consequently enriching the particles' ensemble.
To better analyze the GPS contribution, we propose to compute the prior and posterior GPS errors on the most homogeneous MODIS pixels of the images in terms of land cover, in order to assess the class temperature errors.Since, at the kilometric scale and in the studied region, the landscape is quite heterogeneous, we selected all of the pixels presenting a fraction of the majority class larger than 55%.
Table 4 presents the results obtained in terms of biases, RMSE, and E ff for the seven considered classes.The number of selected pixels is also indicated in the table.The overall statistics show the same features as shown previously, i.e., a small impact of the assimilation window length and the best results obtained for Experiment 1.The biases appear to have been strongly reduced for all classes, as well as the RMSEs.The best improvements are observed for the less frequent class (grasslands) with a reduction of the RMSE of 3.7 K (efficiency index of 67%).The water class is the one for which the assimilation is the less performing even though the efficiency criteria is equal to 30%.This is probably linked to the larger model errors simulated for this class not correctly simulated by the SETHYS model, as previously mentioned.Indeed, a poor model spread was obtained for the prescribed range of variation of the parameters for the water class (not shown here).In other words, the solution space, being outside of the model simulated space, cannot be explored by the GPS and the selected particles remain far from the optimal solution.Nevertheless, on the whole, the results appear quite satisfactory and show the significant contribution of the Meteosat LST data combined with the LSM to estimate relevant subscale temperatures.Concerning the impact of the assimilation window length, our results show better performances for the shortest time period (one day); this was not expected initially.In previous works ( [11,19]), we indeed found that thermal-hydraulic parameters were better calibrated if several days were observed, especially after a rainfall event.In our study case, it appears that radiative parameters were the most sensitive and present a large uncertainty.Therefore, it is comprehensible that a frequent update as provided by a day-time assimilation window, allowing a better fit the observations.

GPS Downscaling Approach Validation at 60 m Scale
Similarly to the previous evaluation at the MODIS scale, a comparison has been performed at 60 m spatial resolution (original resolution of Landsat 7 LST data).The same statistics (RMSE, efficiency, and bias) have been calculated for the same downscaling experiments, over our set of seven Landsat images and by land cover class.It should be noted that at the Landsat 7 scale, each pixel is considered homogeneous.
As an example, Figure 3 presents the results obtained for DOY 240 at 60 m resolution, on the same extracted zone than previously (similarly to Figure 2).We should point out that the GPS downscaled LST image (at 10:15 UTC) and the Meteosat one present a slight delay of 5 min compared to the Landsat 7 image acquired at 10:20 UTC, which is still acceptable given the temporal dynamics of LST.The minimum and maximum LST values observed or estimated by the GPS for the whole image (and not only the extracted zone), are 291.3K and 333 K for Landsat 7, and 297.7 K and 320 K for the GPS.As indicated, the downscaled image does not cover the range of variation of the Landsat one: 21 K compared to 41 K, highlighting the fact that the search for an average temperature per class is probably insufficient at this scale, given the high spatial heterogeneity of the landscape and the intra class variability not accounted by the methodology.However, the large spatial variations are correctly retrieved, even if the large bare soil areas in the southern part of the sub-image appear warmer in the GPS-sharpened image, and the water and some of the winter crop areas appear cooler than in the Landsat product.Some large errors seen on some small agricultural fields give the suspicion of land cover map errors or irrigation practices which were not accounted for in the modeling, the atmospheric forcing data (and the water input) being unique for each land cover type considered.
m spatial resolution (original resolution of Landsat 7 LST data).The same statistics (RMSE, efficiency, and bias) have been calculated for the same downscaling experiments, over our set of seven Landsat images and by land cover class.It should be noted that at the Landsat 7 scale, each pixel is considered homogeneous.
As an example, Figure 3 presents the results obtained for DOY 240 at 60 m resolution, on the same extracted zone than previously (similarly to Figure 2).We should point out that the GPS downscaled LST image (at 10:15 UTC) and the Meteosat one present a slight delay of 5 min compared to the Landsat 7 image acquired at 10:20 UTC, which is still acceptable given the temporal dynamics of LST.The minimum and maximum LST values observed or estimated by the GPS for the whole image (and not only the extracted zone), are 291.3K and 333 K for Landsat 7, and 297.7 K and 320 K for the GPS.As indicated, the downscaled image does not cover the range of variation of the Landsat one: 21 K compared to 41 K, highlighting the fact that the search for an average temperature per class is probably insufficient at this scale, given the high spatial heterogeneity of the landscape and the intra class variability not accounted by the methodology.However, the large spatial variations are correctly retrieved, even if the large bare soil areas in the southern part of the sub-image appear warmer in the GPS-sharpened image, and the water and some of the winter crop areas appear cooler than in the Landsat product.Some large errors seen on some small agricultural fields give the suspicion of land cover map errors or irrigation practices which were not accounted for in the modeling, the atmospheric forcing data (and the water input) being unique for each land cover type considered.Table 5 presents the overall statistics on the whole set of Landsat images for all the downscaling experiments.The best results are obtained for the Experiment 1 (1-day assimilation window, noised forcing data) like for the MODIS comparison.A systematic positive bias ranging between 1.3 and 1.7 K appears for all the downscaling experiments, which was not present in the prior estimations.The Table 5 presents the overall statistics on the whole set of Landsat images for all the downscaling experiments.The best results are obtained for the Experiment 1 (1-day assimilation window, noised forcing data) like for the MODIS comparison.A systematic positive bias ranging between 1.3 and 1.7 K appears for all the downscaling experiments, which was not present in the prior estimations.The posterior RMSEs are larger compared to the prior estimations with differences ranging between 0.4 K and 0.9 K and the efficiencies are all negative.Table 6 presents the statistics calculated for each land cover class, but only for the best scenario case (Experiment 1), since not much difference was obtained among the scenarios.They allow to better understand the performances of the GPS.The table shows that the worse results are obtained for the bare soil and the tree cover classes for which the efficiencies are both negative especially for the bare soil class where it reaches ´55.6%, corresponding to an RMSE of 7 K and a bias of 5.3 K.The rice, water, summer, and winter crops show a slight improvement compared to the prior estimations, with RMSEs slightly decreased (by about 0.5 K) as well as the biases.However, for all of the other classes, the biases when positive for the prior estimations are increased, the negative biases being decreased.

Discussion
To understand the origins of these differences, we have compared the Meteosat and Landsat 7 LST products at the same scales (on the Meteosat grid) and calculated for each dates, the cumulative distribution functions (CDF) of the LST products.An example is given in Figure 4 for the same previous date (DOY 240).The figure shows clearly that Meteosat data presents a warm bias compared to Landsat 7 of about 4 K.This warm bias explains the results shown in Figure 4 where the downscaled LST of bare soils appear warmer compared to Landsat.They also point out that our methodology could be improved in various directions.First, since it relies on the precise knowledge of the land cover mapping, which is more crucial at fine scales, we hope to improve this point in the near future by using new tools and high-resolution data, such as those provided by the Sentinel 2 program.Automatic algorithms, such as the one proposed by [43], could bring significant accuracy in land cover mapping.Another direction of improvement is to relax the land cover map in order to bring corrections regularly to account for changes and for irrigation practices.The assimilation of solar reflectances and microwave data (e.g., radar backscattering coefficients related to soil moisture) should help to correct these issues.The assumption that LST is isothermal for one land cover class within the coarse-resolution pixel, could also be relaxed by working with probability distribution functions instead of single values.Such a framework could be easily implemented in the GPS to represent the intra-class variability and its time evolution.Finally, since the downscaled LST rely strongly on the prior model estimations and parameter range definition, and given the improvements brought by the enlargement of the simulation ensembles (e.g., by noising the input atmospheric variables), more work is needed to adjust these parameters to account for the atmosphere spatial variability.High-resolution reanalysis, such as that provided by the SAFRAN system [44], are also certainly a track to explore in order to improve meteorological input forcing data.

Conclusions
The methodology presented here allows downscaling of Meteosat high-frequency LST time series at higher spatial resolution.The approach is based on the assimilation of low spatial resolution (LSR) observations in a LSM able to simulate the temperature of the various land cover types (components) inside the Meteosat mixed pixel and the composite thermal emission.The method assumes that LST is stationary for each land cover class inside the Meteosat 3 km × 5 km pixel.It requires the availability of a high-resolution land cover map providing the fractions of each land cover type, as well as atmospheric and surface forcing data to simulate prior estimations of the component temperatures.
The method was applied over a 4539 km 2 area of Southeastern France during a seven-month period in 2009 to downscale Meteosat temperatures to 1 km scale as a first step, using MODIS LST for the validation and, at 60 m scale, to match Landsat 7 data.The various downscaling experiments performed demonstrate the benefits of the GPS approach to improve the prior model simulations and the increased performance when short time assimilation period (one day) is prescribed and The same comparisons performed on each of the Landsat images showed that these deviations are not consistent from one day to the next, but also depend on the landscape, and on the land cover.For example, forested areas appear to present larger warmer biases than bare soils.The inconsistent character of these discrepancies, depending on time of day, season, land cover, etc., and the limited dataset did not allow us to conclude about the respective role of directional effects (that are known to be significant at high spatial resolutions ( [42]), atmospheric corrections, or sensor inter-calibration issues.No satisfactory solution for the quantitative evaluation of the downscaled LST product was then possible.Nevertheless, the results obtained, even if they are more qualitative than quantitative, show the potentiality of the downscaling method up to hectometric scales, provided the availability of a robust land cover map.
They also point out that our methodology could be improved in various directions.First, since it relies on the precise knowledge of the land cover mapping, which is more crucial at fine scales, we hope to improve this point in the near future by using new tools and high-resolution data, such as those provided by the Sentinel 2 program.Automatic algorithms, such as the one proposed by [43], could bring significant accuracy in land cover mapping.Another direction of improvement is to relax the land cover map in order to bring corrections regularly to account for changes and for irrigation practices.The assimilation of solar reflectances and microwave data (e.g., radar backscattering coefficients related to soil moisture) should help to correct these issues.The assumption that LST is isothermal for one land cover class within the coarse-resolution pixel, could also be relaxed by working with probability distribution functions instead of single values.Such a framework could be easily implemented in the GPS to represent the intra-class variability and its time evolution.Finally, since the downscaled LST rely strongly on the prior model estimations and parameter range definition, and given the improvements brought by the enlargement of the simulation ensembles (e.g., by noising the input atmospheric variables), more work is needed to adjust these parameters to account for the atmosphere spatial variability.High-resolution reanalysis, such as that provided by the SAFRAN system [44], are also certainly a track to explore in order to improve meteorological input forcing data.

Conclusions
The methodology presented here allows downscaling of Meteosat high-frequency LST time series at higher spatial resolution.The approach is based on the assimilation of low spatial resolution (LSR) observations in a LSM able to simulate the temperature of the various land cover types (components) inside the Meteosat mixed pixel and the composite thermal emission.The method assumes that LST is stationary for each land cover class inside the Meteosat 3 km ˆ5 km pixel.It requires the availability of a high-resolution land cover map providing the fractions of each land cover type, as well as atmospheric and surface forcing data to simulate prior estimations of the component temperatures.
The method was applied over a 4539 km 2 area of Southeastern France during a seven-month period in 2009 to downscale Meteosat temperatures to 1 km scale as a first step, using MODIS LST for the validation and, at 60 m scale, to match Landsat 7 data.The various downscaling experiments performed demonstrate the benefits of the GPS approach to improve the prior model simulations and the increased performance when short time assimilation period (one day) is prescribed and uncertainties in the atmospheric forcing are accounted for.The results obtained at 1 km scale, show that the downscaled temperatures can be estimated with an uncertainty better than 2.7 K which is an improvement compared to the prior model values of 4 K.Further comparison to smaller scales could not be quantitatively evaluated because of the lack of unbiased data for the validation.Finally, the crucial knowledge of the land cover map was highlighted and improvement tracks, like the coupled assimilation of auxiliary remote sensing products, were suggested.
This study aims to demonstrate the relevance of data assimilation approaches, and especially GPS, to perform LST downscaling.The method requires the implementation of a LSM at finer scales and the various input datasets needed to run the model.Given that, it allows dynamic monitoring of the subpixel LSTs at high spatio-temporal resolutions.As a result, the method not only provides downscaled maps at the time of observation, but at the model time step (in our case, each half-hour) and, therefore, even in cloudy conditions.These results also show the potentialities of the GPS for the assimilation of remote sensing data over mixed pixels.If the dynamic model is able to simulate the composite signal, the GPS can be implemented efficiently to optimize subpixel model parameters.The added value compared to other DA approaches is the easy development and parallel implementation of the downscaling process that could present an interesting framework for an automatic real-time implementation.

Figure 1 .
Figure 1.Crau-Camargue region and experimental sites equipped with micrometeorological measurements (upper panel).The white box delineates the region under study corresponding to the Landsat 7 frame.The land cover map of the study area is plotted in the lower panel.The red box indicates the region highlighted in the validation study.

Figure 1 .
Figure 1.Crau-Camargue region and experimental sites equipped with micrometeorological measurements (upper panel).The white box delineates the region under study corresponding to the Landsat 7 frame.The land cover map of the study area is plotted in the lower panel.The red box indicates the region highlighted in the validation study.

Figure 2 .
Figure 2. GPS downscaling results at 1km resolution over the study area and for the DOY 240: (a) land cover dominant class; (b) Meteosat LST; (c) GPS LST; (d) MODIS LST; and (e) map of the differences between GPS and MODIS.Both Meteosat and MODIS LST were acquired at 11:00 UTC and GPS results are obtained at the same time.

Figure 2 .
Figure 2. GPS downscaling results at 1km resolution over the study area and for the DOY 240: (a) land cover dominant class; (b) Meteosat LST; (c) GPS LST; (d) MODIS LST; and (e) map of the differences between GPS and MODIS.Both Meteosat and MODIS LST were acquired at 11:00 UTC and GPS results are obtained at the same time.

Figure 3 .
Figure 3. GPS downscaling results at 60 m resolution over a small part of the studied area for the DOY 240: (a) land cover dominant class; (b) Meteosat LST; (c) GPS LST; (d) Landsat 7 LST; and (e) map of the differences between GPS and Landsat 7. Meteosat data was acquired at 10:15 UTC and Landsat 7 at 10:20 UTC.

Figure 3 .
Figure 3. GPS downscaling results at 60 m resolution over a small part of the studied area for the DOY 240: (a) land cover dominant class; (b) Meteosat LST; (c) GPS LST; (d) Landsat 7 LST; and (e) map of the differences between GPS and Landsat 7. Meteosat data was acquired at 10:15 UTC and Landsat 7 at 10:20 UTC.

Figure 4 .
Figure 4. Cumulative distribution function (CDF) calculated for the three images acquired at about the same time by Meteosat (Blue curve) and Landsat 7 upscaled at the same resolution (green curve) and Landsat 7 at full resolution (30 m, pink curve).

Figure 4 .
Figure 4. Cumulative distribution function (CDF) calculated for the three images acquired at about the same time by Meteosat (Blue curve) and Landsat 7 upscaled at the same resolution (green curve) and Landsat 7 at full resolution (30 m, pink curve).

Table 1 .
Description of the selected parameters for the considered land cover types.

Table 2 .
Experimental settings for the considered experiments.

Table 3 .
GPS downscaling results averaged over the whole region and for the 23 considered validation dates, compared to SETHYS prior estimations.

Table 3 .
GPS downscaling results averaged over the whole region and for the 23 considered validation dates, compared to SETHYS prior estimations.

Table 4 .
GPS downscaling results by class averaged over MODIS homogeneous pixels and the 23 considered validation dates, compared to SETHYS prior estimations.

Table 5 .
GPS results averaged over the whole region and for the seven considered validation dates, compared to SETHYS prior estimations.

Table 6 .
GPS downscaling results by class averaged over the seven considered validation dates at 60 m resolution, compared to SETHYS prior estimations.