Kalman filter physical retrieval of surface emissivity and temperature from geostationary infrared radiances

The high temporal resolution of data acquisition by geostationary satellites and their capability to resolve the diurnal cycle allows for the retrieval of a valuable source of information about geophysical parameters. In this paper, we implement a Kalman filter approach to apply temporal constraints on the retrieval of surface emissivity and temperature from radiance measurements made from geostationary platforms. Although we consider a case study in which we apply a strictly temporal constraint alone, the methodology will be presented in its general four-dimensional, i.e., space-time, setting. The case study we consider is the retrieval of emissivity and surface temperature from SEVIRI (Spinning Enhanced Visible and Infrared Imager) observations over a target area encompassing the Iberian Peninsula and northwestern Africa. The retrievals are then compared with in situ data and other similar satellite products. Our findings show that the Kalman filter strategy can simultaneously retrieve surface emissivity and temperature with an accuracy of ± 0.005 and±0.2 K, respectively.


Introduction
Infrared instrumentation on geostationary satellites is now rapidly approaching the spectral quality and accuracy of modern sensors on board polar platforms.Currently at the core of European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) geostationary meteorological programme is the Meteosat (meteorological satellite) Second Generation (MSG).EUMETSAT is preparing for Meteosat Third Generation (MTG), which will carry the Flexible Combined Imager (FCI) with a spatial resolution of 1-2 km at the sub-satellite point and 16 channels (8 in the thermal band), and an infrared sounder (IRS) that will be able to provide unprecedented information on horizontally, vertically, and temporally (four-dimensional; 4-D) resolved water vapor and temperature structures of the atmosphere.The IRS will have a hyperspectral resolution of 0.625 cm −1 wave numbers, will take measurements in two bands, the long-wave infrared (LWIR) (14.3-8.3 µm) and the mid-wave infrared (MWIR) (6.25-4.6 µm), and with a spatial resolution of 4 km and a repeat cycle of 60 min.
Because geostationary satellites are capable of resolving the diurnal cycle, and hence providing time-resolved sequences or times series of observations, they are a source of information which can suitably constrain the derivation of geophysical parameters.In this paper, we implement a Kalman filter (KF) approach for applying temporal constraints on the retrieval of surface emissivity and temperature from radiance measurements made from MSG SEVIRI (Spinning Enhanced Visible and Infrared Imager).The study has been performed also in view of future applications to the MTG mission.This mission should improve sounding density, quality and accuracy of surface and atmospheric parameters.
The Kalman filter (Kalman, 1960;Kalman and Bucy, 1961) has received widespread attention in the context of numerical weather prediction (NWP) and in the broad research area of data assimilation (e.g., Lorenc, 1986;Evensen, 1994;Talagrand, 1997;Nychka and Anderson, 2010).Specific Published by Copernicus Publications on behalf of the European Geosciences Union.
applications to atmospheric chemistry with satellite data have been described by, e.g., Khattatov et al. (1999); Lamarque et al. (1999); Levelt et al. (1998).For a detailed review and tutorial of the theoretical background of KF the reader is referred to, e.g., Wikle and Berliner (2007); Nychka and Anderson (2010).
The present paper addresses the capability of the Kalman filter to convey temporal constraint in the retrieval of surface parameters though time series of geostationary satellite data.The fact that time continuity of the observations brings much information about atmospheric processes is normally evidenced by the pronounced dynamical correlation, which in many instances can be modeled with Markov chains or Markov stochastic processes (e.g., Serio, 1992;Cuomo et al., 1994;Serio, 1994, and references therein).
To exploit the temporal information, we focus on implementing the KF so that it is capable of incorporating dynamical correlation within the retrieval process without making use of a full dynamical NWP system.We aim at developing a retrieval strategy for surface emissivity and temperature which allows us insight into understanding how we can better exploit satellite data per se.In other words, the analysis is conducted within a context which envisages an almost entirely data-driven system.In this respect, we clarify that although we try to exploit tools such as the KF which are generally used in a data assimilation context, we aim at addressing a retrieval problem limited to surface parameters.We do not want to solve an assimilation problem according to the common method (e.g., see Nychka and Anderson, 2010) of combining a NWP model with observations.
In view of future MTG applications, the KF methodology is presented in a general context which applies to both spatial and temporal constraints.However, it will be exemplified for a case study in which we consider a strictly temporal constraint alone.As said, this is the problem of surface temperature (T s ) and emissivity ( ) separation, that is, the simultaneous retrieval of (T s , ) from SEVIRI infrared channels.Toward this objective, a case study has been defined which includes a specific target area characterized by a large variety of surface features.
The problem of retrieving surface emissivity and temperature from satellite data has long been studied.A recent review of the subject has been provided by Li et al. (2013).According to this review, our KF approach from a geostationary platform is novel.A similarity could be found with the scheme developed by Li et al. (2011).However, while in Li et al. (2011) the observations are accumulated for a prescribed time slot (normally six hours), we pursue a genuine dynamical strategy which exploits the sequential approach of the Kalman filter.This results in an algorithm which does not need to increase the dimensionality of the data space because of time accumulation, while preserving the highest time resolution prescribed by the repeat time of the geostationary instrumentation (15 min for SEVIRI).
To check the quality and performance of our approach, KF retrieval results will be compared with in situ data, spacetime collocated ECMWF (European Centre for Medium-Range Weather Forecasts) analysis and other similar satellite products, such as those from the AVHRR (Advanced Very High Resolution Radiometer).
The study is organized as follows.Section 2 will present the data used in the analysis.This section will also provide some details about the forward model we have developed for SEVIRI.Section 3 will deal with the retrieval methodology, whereas Sect. 4 will exemplify the application of the methodology to a SEVIRI case study.Finally, conclusions will be made in Sect. 5.

Data and forward modeling
In this paper, the KF methodology will be applied for the retrieval of surface emissivities and surface temperature from SEVIRI infrared channels in the atmospheric window over a target area, covering a geographic region with very different surface features: sea water, arid, vegetated and cultivated land, and urban areas.
The SEVIRI imager on board Meteosat-9 allows for a complete image scan (full Earth scan) once every 15 min period with a spatial resolution of 3 km for 12 channels (8 in the thermal band), over the full disk covering Europe, Africa and part of South America.
SEVIRI infrared channels range from 3.9 µm to 12 µm.Their conventional definition in terms of channel number is given in Table 1, whereas their spectral response is shown in Fig. 1.The figure also provides a comparison with a typical IASI (Infrared Atmospheric Sounder Interferometer) spectrum at a spectral sampling of 0.25 cm −1 .
IASI has been developed in France by the Centre National d'Etudes Spatiales (CNES) and is on board the Metop (Meteorological Operational Satellite) platform, a series of three satellites belonging to the EUMETSAT European Polar System (EPS).The instrument has a spectral coverage extending from 645 to 2760 cm −1 , which, with a sampling interval σ = 0.25 cm −1 , gives 8461 data points or channels for each single spectrum.Data samples are taken at intervals of 25 km along and across track, each sample having a minimum diameter of about 12 km.Further details on IASI and its mission objectives can be found in Hilton et al. (2012).Atmospheric parameters (temperature, water vapor and ozone profiles) derived from IASI spectral radiances will be used in this paper to assess the sensitivity of SEVIRI atmospheric window infrared channels to the atmospheric state vector.
As stated before, for the purpose of this study, SEVIRI Meteosat-9 high rate level 1.5 image data and IASI (level 1C) observations have been collected for the target area shown in Fig. 2 for the full month of July 2010.The area is covered with 392 088 Meteosat-9 pixels and includes Spain, Portugal,  part of the northwestern Africa, and the western part of the Mediterranean Basin.
To check the performance of the scheme, we have also selected three smaller areas (also shown in Fig. 2 with red boxes) in Spain, the Sahara desert, and the Mediterranean Basin, which have a size of 0.5 × 0.5 degrees and each correspond to one box of the ECMWF analysis grid mesh (e.g., see Fig. 3).For the Spanish location the area includes 187 SEVIRI pixels, 219 for the Sahara desert and 178 for the Mediterranean Basin.The land coverage for the small target area close to Seville is a mosaic of cultivated areas, with green grass, foliage, bare soil and urban areas.For this type of coverage we expect an emissivity at atmospheric window well above 0.90.The small Sahara desert area is just a desert sand homogeneous flat area, with no vegetation.In this case we know that emissivity is dominated by quartz particles, which yield a characteristic fingerprint at 8.6 µm (reststrahlen doublet of quartz).This strong signature is in the middle of the SEVIRI channel at 8.7 µm, and, therefore, the retrieved emissivity at this channel has to show the quartz fingerprint.
observations, an introduction.J. Me-9, 1997.Theory: Methods for Data Fitting ation.Elsevier, New York, 1987. esen, F., and Kabsch, E.   For the whole target area shown in Fig. 2, we have also acquired ancillary information for the characterization of the thermodynamical atmospheric state.This information is provided by ECMWF analysis products for the surface temperature, T s and the atmospheric profiles of temperature, water vapor and ozone (T , Q, O) at the canonical hours 00:00, 06:00, 12:00 and 18:00 UTC.ECMWF model data are provided on 0.5 × 0.5 degree grid.In each ECMWF grid box there are on average ≈ 200 SEVIRI pixels, for which we assume that the atmospheric state vector is the time collocated ECMWF analysis (e.g., see Fig. 3).
Within the inverse scheme, an important issue concerns a priori information to constrain the retrieval of emissivity.To this end, we have used the University of Wisconsin Baseline Fit Global Infrared Land Surface Emissivity Database (UW/BFEMIS database, e.g., http://cimss.ssec.wisc.edu/iremis/)(Seemann et al., 2008;Borbas and Ruston, 2010).The UW/BFEMIS database is available for years 2003 to 2012, globally, with 0.05 degree spatial resolution.Details of how to transform UW/BFEMIS database emissivity to SE-VIRI channel emissivity can be found in Serio et al. (2013); Masiello et al. (2013).
For the purpose of comparison, we have used also NOAA (National Ocean and Atmosphere Administration) Optimum Interpolation 1/4 degree daily Sea-Surface Temperature (OISST) analyses for the month of July 2010.The analysis, which is a product of the processing of AMSR (Advanced Microwave Scanning Radiometer) and AVHRR, will be compared to that obtained by SEVIRI for sea surface.The    analysis will be referred to as AMSR+AVHRR OISST in the remainder of this paper.The AMSR+AVHRR OISST analysis has been downloaded from the website ftp://eclipse.ncdc.noaa.gov/pub/OI-daily-v2/NetCDF/2010/AVHRR-AMSR/.
Finally, we have also used data collected at the Evora ground site (38.55 • N, 8.01 • W) located in southern Portugal and maintained by the EUMETSAT Satellite Applications Facility on Land Surface Analysis (LSA SAF) team.The area surrounding the site is dominated by Quercus woodland plains and is fairly homogeneous at the SEVIRI spatial scales (Dash et al., 2004).The ground station is equipped with a suite of radiometers (9.6-11.5 µm range) providing temperatures of tree canopies and of ground in the sun and in the shade.These are combined to provide a composite ground temperature representative of SEVIRI pixels, considering that the fractional area coverage of canopies is 0.32 (Trigo et al., 2008).It is worth mentioning that the ground and the treetop canopy present contrasting temperatures particularly during daytime, when differences can easily reach 15 K.As a consequence, the composite ground temperatures are fairly sensitive to the fraction of trees being considered.For this purpose, the area surrounding the Evora station was carefully characterized using very high resolution IKONOS satellite images (Kabsch et al., 2008;Trigo et al., 2008).In addition, our Kalman Filter retrievals for the pixels closer to Evora are also compared with the operational land-surface temperature product provided by the LSA SAF (Freitas et al., 2010).

Forward models: σ -IASI and σ -SEVIRI
Forward calculations for the SEVIRI channels 2-8 (see Table 1) are obtained by the σ -SEVIRI code that we have developed specifically for this study.
We do not consider the SEVIRI channel at 3.9 µm since during daytime it is contaminated by reflected solar radiation and affected by non-local thermodynamic equilibrium (non-LTE) effects.Furthermore, the CO 2 line mixing at 4.3 µm CO 2 band head is poorly modeled in state-of-the-art radiative transfer and can add potentially large bias.
Regarding channels 2 to 8, the forward model, σ -SEVIRI has been derived from σ -IASI (Amato et al., 2002) which is a monochromatic radiative transfer designed for the fast computation of spectral radiance and its derivatives (Jacobian) with respect to a given set of geophysical parameters.
The form of the radiative transfer equation, which σ -IASI and hence σ -SEVIRI consider in its numerical scheme, has been recently reviewed and presented in Masiello and Serio (2013), to which the interested reader is referred.The model also takes into account the radiance term, which is the radiation reflected from the surface back to the satellite.Both Lambertian and specular reflections can be modeled.
To accomplish the radiative transfer calculation σ -IASI uses a lookup table for the optical depth; this table was developed from one of the most popular line-by-line forward models, Line-By-Line Radiative Transfer Model (LBLRTM) (Clough et al., 2005).
The model σ -SEVIRI is itself based on a lookup table, which is obtained by a proper down-sampling of the lookup table for σ -IASI.For this reason we need to give some details about σ -IASI in order to describe how σ -SEVIRI works.
The σ -IASI model (Amato et al., 2002) parameterizes the monochromatic optical depth with a second-order polynomial.At a given pressure-layer and wave number σ (in cm −1 units), the optical depth for the a generic ith molecule is computed according to where T is the temperature, q i the molecule concentration and c σ,j,i with j = 0, 1, 2, are fitted coefficients, which are actually stored in the optical depth lookup table.
For water vapor, unlike other gases, in order to take into account effects depending on the gas concentration, such as self-broadening, a bi-dimensional lookup table created by Masiello and Serio (2003) is used.Thus, for water vapor, identified with i = 1, the optical depth is calculated according to The subscript σ indicates the monochromatic quantities.In the case of hyperspectral instruments, such as IASI, the monochromatic optical depths are computed and parameterized at the spectral sampling interval of 10 −4 cm −1 .This spectral sampling is much too fine for a band instrument such as SEVIRI.In the case of SEVIRI, the spectral Atmos.Meas.Tech., 6, 3613-3634, 2013 www.atmos-meas-tech.net/6/3613/2013/sampling can be averaged and sampled at a rate of 10 −1 cm −1 without sacrificing accuracy.Also in this case the optical depth can be parameterized with a low-order polynomial, and its coefficients are obtained as explained below.
For each species, i, we can define an equivalent optical depth which can be parameterized with respect to temperature in the same way we do for monochromatic quantities (Eqs. 1 and 2).Considering the larger channel bandwidths of the SEVIRI measurements, averaging is applied over the spectral wave-number band of each channel.This averaging is identified by the angular brackets • .The equivalent optical depth is where the equivalent coefficients c σ ,j,i with j = 0, 1, 2, are obtained by fitting the layer transmittance averaged over the coarse sampling of 10 −1 cm −1 , Because of this down-sampling σ -SEVIRI, which is based on the coarse-mesh lookup table, runs ≈ 1000 times faster than σ -IASI.As the parent code, σ -IASI, σ -SEVIRI can compute the analytical Jacobian derivative for a large set of surface and atmospheric parameters: , T s and (T , Q, O).

The retrieval framework
Before showing the retrieval problem for the pair of surface parameters (T s , ), we briefly review the concept of the optimal estimation in the general context of data assimilation (Lorenc, 1986;Talagrand, 1997;Wikle and Berliner, 2007;Nychka and Anderson, 2010;Rodgers, 2000), which allows us to describe the retrieval methodology in its general spatiotemporal framework and also to put in evidence its commonalities with the KF methodology.
For the benefit of the reader, we will try to stay as close as possible to the notation used in Rodgers (2000), therefore the symbol and subscript ε will be used to denote the observational covariance matrix and hence the noise term affecting the spectral radiance.For emissivity we will use the symbol , which should not be confused with ε.

Static a priori background
To simplify the exposition let us assume that the times are indexed by integers, t = 1, 2, . .., although handling unevenly spaced times does not add any fundamental difficulty.The derivation of the thermodynamical state of the atmosphere, at a given time t, given a set of independent observations of the spectral radiance, R t (σ ), is well established when each time t is considered independent from past and future measures.Let R t be the radiance vector with m the number of spectral radiances, and where the superscript T stands for transpose.Under the assumption of multivariate normality the retrieval problem can be seen as one of variational analysis in which a suitable estimation of the state vector is obtained by minimizing the form (see e.g., Courtier, 1997;Talagrand, 1997;Tarantola, 1987;Carissimo et al., 2005): where F is the forward model function; v is the atmospheric state vector, of size n; v a is the atmospheric background state vector, of size n; S ε is the observational covariance matrix, of size m × m; and S a is the background covariance matrix, of size n × n.Equation ( 6) is commonly linearized and a Gauss Newton iterative method is used to solve the quadratic form where It should be stressed that, formally, the state vector, v can be thought of as a 3-D geophysical field, and not necessarily of a vector in one dimension (altitude coordinate).
In the context of data assimilation, x a is normally the forecast at time t, and S a is the forecast error covariance matrix.The estimation, x, is referred to as the analysis.

The Kalman filter
The Kalman filter was first developed by Kalman (1960) and Kalman and Bucy (1961) in an engineering context and as a linear filter.Its derivation from the Bayes formalism has been shown by many authors (e.g., see the review by Wikle and Berliner, 2007).With our notation, the formal filter can be summarized with the two equations below which are often referred to as the observation equation (or data model) and the state equation (or dynamic model or system model), respectively.Here M is a linear operator and the noise model term, η t has covariance, S η .The remaining parameters appearing in Eq. ( 9) have the same meaning as those introduced in Sect.3.1.KF is intrinsically linear, therefore the observation equation has to be linearized in order to write down the optimal estimation for the state vector.With the same notation we have used until now, we have the linear KF form where we use the notation K t for the Jacobian to stress that it depends on time, t.
It should be noted that we assume that both the noise terms ε t and η t are independent of the state vector.

The KF update step or analysis
Under the same assumption of multivariate normal statistics as that used in Sect.3.1, we have that the optimal KF estimate, xt at time t is given by (e.g., Wikle and Berliner, 2007) We see that the optimal KF estimate for xt is formally equivalent to that obtained by the variational or optimal estimation approach in Sect.3.1.We recall, once again, that in the context of data assimilation, x a is normally the forecast at time t, and S a is the error forecast covariance matrix.The estimation, xt , is referred to as the analysis at time t, which has covariance matrix given by Ŝt .
One important aspect of the formal solution is that the analysis update depends only on the data at time t and not on that at previous times.This property is referred to as the Markov property.In fact, the formal solution for the analysis does not depend on the dynamical system directly.We can see that the expression in Eq. ( 11) does not contain the operator M.
The above property is also referred to as the regularization property of KF.New data comes in at t and the KF updated state estimate is the minimizer of the quadratic form or cost function, S: However, an important distinction regarding data assimilation is that S a is potentially generated from the process and not from an external spatial model.In fact S a is iterated with the process, as will become clear in examining the forecast step for the linear KF.It is important here to stress that the minimization of the form (12) needs an iterative approach because of the nonlinearity of the forward model and a criterion to stop iterations.We use the usual χ 2 criterion.In fact, under linearity, the value of twice the quadratic S (Eq.12) at the minimum is distributed as a χ 2 variable with m degrees of freedom (Tarantola, 1987).A χ 2 threshold, χ 2 th , at three sigma confidence intervals, can be then obtained according to χ 2 th = m + 3 √ 2m.Therefore, the iterative procedure is stopped when

The KF forecast step
In our notation, x = v − v o and xa = va − v o , so that the formal KF estimate for the state vector is For the forecast step the KF assumes that the process evolves in a linear way, according to the operator M; therefore, we can obtain an estimate of the forecast at time t + 1, standing at time t, through the linear transform where the superscript f stands for forecast.The forecast has uncertainty given by where S η is the covariance matrix of the noise term η t (see Eq. 10).
As soon as new data comes in at time t + 1, the forecast becomes the background, and we are ready to obtain the new analysis, vt+1 .An important concept to draw from this sequential updating is that spatial information about the distribution of v t can be generated from the dynamics of the process.In fact, analyzing the forecast covariance matrix (Eq.16), it is seen that it is based on the previous forecast covariance matrix and also inherits the dynamical relationship from the previous time.Thus, in the situation of assimilation for a space-time process, the spatial covariance for inference is built up sequentially based on past updates with observations and propagating the posterior forward in time as a forecast distribution.We stress that this spatial information is the difference or error between the conditional mean and the true field and is not the covariance of the process itself.
However, the goodness of this spatial information mostly relies on the quality of the physics we model with the operator M. Typically, the forecast step is completed by a deterministic, physically based model.In this case, the spatial information has value.However, in a case in which we want Atmos.Meas.Tech., 6, 3613-3634, 2013 www.atmos-meas-tech.net/6/3613/2013/ the problem driven from the data, the model can be very simplistic and inherently inadequate to describe the real-world process.In this case, spatial information has to be provided externally through a proper definition of S a .

A formulation of the emissivity/temperature retrieval with KF
As stated at the beginning of Sect.3, the general KF formalism has been described and presented in a full 4-D setting.
In this section we will deal with an application to the (T s , ) problem where we apply a strictly temporal only method.
To begin with, we introduce a transform for the emissivity, which allows us to constrain the retrieval to the physical emissivity range of 0-1.Letting be the emissivity at any of the channels, we consider the logit transform which has inverse The transform maps 0-1 into the interval [−∞, +∞] and vice versa.Therefore if we work with the variable e, retrieval positiveness for is ensured.
In order to work with the parameter e we have to properly transform the Jacobian.It easily follows from Eq. ( 18) that where R is the radiance at a generic channel.
If we linearize the forward model, at time t, with respect to e and T s , we obtain with δe = e − e o of dimension m × 1, δT st = T st − T sto .The matrix A t is the emissivity Jacobian, a diagonal matrix of size m × m, and B t is the surface temperature Jacobian, a vector of dimension m × 1.We have that the size of the observation vector, y t is m × 1, the dimension of the Jacobian , and that the state vector, has dimension (m + 1) × 1.As regards the state or model equation for emissivity, an evolution equation is straightforward if we consider the high repeat rate of SEVIRI observations (15 min).This leads us to assume that the evolution of has a low variability on a time scale of few hours.This is particularly true for emissivity, but much less for temperature over land, which is strongly influenced from the daily cycle.For sea surface the assumption of a low time variability on time scales of several hours is good both for T s and .
With this in mind, let v = (e 1 , . .., e m , T s ) T be the emissivity-temperature vector, a suitable dynamical equation is then a simple persistence where, according to our notation (see Sect. 3.2), η t is a noise term with covariance, S η , and M is the identity propagation operator.
We know that the persistence model of Eq. ( 23) is not physically correct since it cannot reproduce the strong daily cyclic behavior of T s expected in clear sky for land surface (Gottsche and Olesen, 2009;Menglin and Dickinson, 1999;Menglin, 2000).It could be a fair model for sea surface, where thermal inertia of water strongly damps the effect of the solar cycle; however, it cannot represent a good model for land surface.
Nevertheless, it has to be stressed that within the context of the Kalman filter methodology we can accommodate our knowledge about the adequacy of the model (Wikle and Berliner, 2007).In practice, provided that the parameters are strongly constrained by the data, the precise form of the evolutionary equation is not important for the estimation problem as long as the error covariance appropriately reflects the uncertainty of the current state estimate.To this end, an important role is played by the stochastic noise covariance, S η .By properly tuning the stochastic noise covariance, we can have a retrieval which is either dominated by the data (S η → +∞, model inadequate), or the state model (S η → 0, model adequate).
SEVIRI atmospheric window channels are strongly dominated by T s .This is exemplified in Fig. 4, which shows a simulation of the daily evolution of T s for a desert site and the corresponding radiance signal at channel 7 (12 µm).The simulation has been obtained using the daily cycle model developed by Gottsche and Olesen (2009) on the basis of in situ observations made at a station in the Namib desert.The model fits the data with an accuracy of ≈ 1-2 K, therefore the T s evolution shown in Fig. 4 reflects a realistic situation.
The corresponding radiance has been obtained through σ -SEVIRI.The state vector needed for the computation of the radiance has been obtained from the ECMWF analysis for a desert site.
Another way to assess the strong dependence of the SE-VIRI atmospheric window channels on temperature is to compute the index D between two consecutive times, t and t + 1, defined according to where σ j denotes the wave number of a generic SEVIRI channel.This index is the ratio of the derivative Jacobian of the surface temperature to the increment of the radiance due to the variation of the surface temperature within the time interval (t, t + 1).Because of the meaning of the Jacobian, this index has to be close to 1 in case the channel strongly depends on T s .Note that the second factor in Eq. ( 24) is the finite-difference-based calculation of the inverse of the Jacobian itself.For the case shown in Fig. 4a, b, Fig. 4c shows the index D for the three SEVIRI atmospheric channels.From this figure it is immediately seen that the radiance time-behavior is completely dominated by the time-evolution of T s .This is a helpful situation because, at least for temperature, we can design a Kalman filter which is strongly driven by the data.
To this end, we first clarify how we build up S η and S a on the basis of the related matrices for emissivity and surface temperature.
We do not consider correlation between emissivity and surface temperature; therefore, and where S ηe is the covariance matrix of the emissivity stochastic term; S ηT s is the variance (scalar) of the surfacetemperature stochastic term; S e is the initial background covariance matrix of the emissivity vector; and S T s is the initial background variance (scalar) of the surface-temperature parameter.
At this point we have defined all the components of our (T s , ) problem which are needed to run the Kalman filter.The flow of operations is here summarized for the benefit of the reader.First, obtain the analysis update through Eq. ( 14); second, compute the forecast with Eq. ( 15); third, find the forecast covariance matrix through Eq. ( 16); and, finally, define the forecast to be the new background (Eq.17) and return to Eq. ( 14) for a new cycle.
Further details of how we build up the above ingredients are given below.To begin with, S e is derived from the UW/BFEMIS database (see Sect. 2).Its definition and calculation is space-time localized.For a given month and SE-VIRI pixel location, UW/BFEMIS yields ten different samples of the emissivity vector from ten different years.The UW/BFEMIS emissivity-vector samples undergo the logit transform (see Eq. 18) and are used to compute the covariance matrix S e .It could be argued that S e built up on ten samples implies an unrealistically large statistical uncertainty for the covariances.This is true and reflects the present stage of our knowledge about surface emissivity.Because of this uncertainty we will be forced to apply somewhat ad hoc further manipulations of this covariance matrix to get realistic retrievals.
An example of S e for the set of seven SEVIRI channels (2 to 8 in Table 1), for the month of July and for a SEVIRI pixel corresponding to a site in the Sahara desert is shown in Table 2. S e shown in Table 2 makes reference to the emissivity vector ordered from longest to shortest wave number.The element (5,5) corresponds to the channel at 8.7 µm, which is in the middle of the quartz reststrahlen band and hence is characterized by the strongest variability.
The covariance matrix S ηe is derived from S e with a scaling procedure.This is justified because of the need to scale down S ηe in order to correctly take into account the expected variation of emissivity on a time scale comparable to the SE-VIRI repeat time of 15 min.However, this is a rather ad hoc inflation/deflation procedure, which is performed on the basis of trial and error until we yield realistic retrievals.
The covariance matrix S e is scaled according to the following procedure.Let S e (i, j ); i, j = 1, . .., m = 7 be the elements of S e .The correlation matrix C e is defined according to and the matrix S e is scaled according to where S (s) e is the matrix scaled by the vector of scaling factors, f .The scaling operation above preserves the correlation structure and in practice we consider a constant scaling factor, f 2 , that does not change along the diagonal.We assume S ηe = S (s)

Atmos
e .As already mentioned, the appropriate value of f has to be tuned in simulation.After extensive simulations (Serio et al., 2013), we have found that f = 10 is appropriate for this case study.
As far as T s is concerned, based on the evidence of Fig. 4, we want to stay closer to the data than to the model.We have that a variance of 1 K 2 for the initial background and stochastic terms S T s and S ηT s respectively, provides a balanced retrieval.In other words, at least for land surface, S ηT s does not need to be downscaled with respect to S T s .
This can be seen in Fig. 5, where we show the results of a retrieval exercise obtained in simulation for the case of desert site (see Serio et al., 2013, for full details).The case shown uses a persistence model for the state equation of both emissivity and skin temperature.For emissivity this is correct, since the simulation assumes a constant emissivity at each SEVIRI channel.Conversely it is not correct for the surface temperature, whose true value follows the daily cycle shown in Fig. 5.
The example shown in Fig. 5 and the error analysis in Fig. 6 allows us to illustrate the property of the KF to accommodate the knowledge of the adequacy of the deterministic model.As said before this is obtained by properly tuning the stochastic term.In the example shown in Fig. 5, the stochastic variance for T s has been set to 1 K 2 .Because of this choice, we correctly follow the data and retrieve the true value of the surface temperature within the accuracy determined by the a posteriori covariance matrix, that is ≈ 0.2 • C. The same conclusion holds for emissivity, which is retrieved within an accuracy of ≈ 0.005.For the stochastic variance, we can also prescribe a value just equal to zero, which will result in a retrieval highly dominated by the model.The results are presented in Fig. 7.We see that after some iterations, the retrieval just follows the (inadequate) persistence model.The initialization point for skin temperature in both exercises is the true temperature minus 4 • C. Note that we need to specify only the initialization point at t = 0, after that KF yields the retrieval on the basis of the data points and model alone.The results shown in Figs. 5 to 7 justify the use of a simplistic persistence model for T s because the data -that is, observations -strongly constrain the phenomenon under investigation.We can get to the same conclusion by observing that, even by prescribing a stochastic variance for T s equal to 1 K 2 , we obtain a precision for the final estimate of ≈ 0.2 • C.This finding implies that the a priori information for T s has little to no impact on the final estimate, which is therefore largely dominated by the data.
The retrieval exercise shown in Figs. 5 to 7 allows us to address another important issue: how the time constraint imposed with the persistence model improves the retrieval in comparison to a scheme where this constraint is not imposed at all.If we try to solve the (T s , ) retrieval problem within the usual context of least squares estimation with a static background, we get a retrieval covariance matrix with a strong anticorrelation between T s and .This anticorrelation is due to the relationship of these two parameters within the radiative transfer equation and makes their effective separation unpractical.To exemplify this effect, we have run the same retrieval exercise shown in Fig. 5, but now with M = 0.In this way, the retrieval does not evolve through the state equation and the surface parameters are estimated on the basis of a static background.The error analysis for this exercise is shown in Fig. 8.The anticorrelation effect is soon evident.We find that T s is biased significantly low and significantly high.Comparing Fig. 6 with Fig. 8, we can clearly identify the impact of temporal information propagation from the KF versus the case without this propagation.The comparison confirms the merit of the KF application for this problem.
It is also noteworthy that for land we have empirical state models, which can reproduce the surface temperature daily cycle with high accuracy (Gottsche and Olesen, 2009;Menglin and Dickinson, 1999;Menglin, 2000).Therefore the question may be posed as to whether or not we can improve the results by using a more adequate model for T s .This exercise has been performed in Serio et al. (2013), where the temperature daily cycle was modeled with a second-order autoregressive process.However, no improvement was found with respect to a simple persistence model.The fact is that the daily cycle is reproduced in its very fine details by the data, as it is possible to see, e.g., from Fig. 4. Therefore, for the particular case of retrieving (T s , ), there is no essential need to include the daily cycle information though an external model.However, it has also to be stressed that this may not be the case, e.g., for atmospheric parameters where the satellite infrared observations are less adequate and need to be assimilated in a system with an accurate dynamical model.
Finally, we stress that for sea surfaces a simple persistence model is accurate also for the skin temperature, and therefore S ηT s needs to be downscaled with respect to S T s .We use S T s = 1 K 2 and obtain S ηT s again by scaling with a factor f = 10, that is S ηT s = 0.01 K 2 .
For the sea-surface emissivity covariance we use Masuda's model (Masuda et al., 1988).For any single SEVIRI  pixel field of view angle, we generate the emissivity vector for wind speed in the range 0-15 m s −1 and with a step of 1.5 m s −1 .In this way we have 11 emissivity vectors, which are used to define background vector and covariance.Again, the resulting covariance is downscaled by a factor f = 10.
The validity of the persistence model for sea surface has been checked directly on the basis of real observations, because for sea surface the ECMWF analysis is credited with an accuracy within ±1 K. Figure 9 exhibits the results for the sea target area shown in Fig. 2 and for 31 July 2010.We see that a stochastic variance term below 0.25 K 2 tends to have a better agreement with the ECMWF model, which leads us to conclude that for sea surface a persistence model is effective not only for emissivity, but also for T s .
In passing, we also note from Fig. 9 that the skin temperature reaches a maximum around 15:00 UTC.The maximum around 15:00 UTC is in agreement with Gentemann et al. (2003), who showed that during the daytime, solar heating may lead to the formation of a near-surface diurnal warm layer, particularly in regions with low wind speeds.Analysis   of TMI (the Tropical Rainfall Measuring Mission (TRMM) Microwave Imager) and AVHRR skin temperature have revealed that the onset of warming begins as early as 08:00 and peaks near 15:00 with a magnitude of 2.8 • C during favorable conditions.

Sensitivity to the atmospheric state vector
For the problem of (T s , ) retrieval, we consider SEVIRI atmospheric window channels alone, namely channels 4, 6, and 7 (see Table 1).However, in practice atmospheric window channels can have a contribution from the atmospheric parameters (T , Q, O) which give the major emission contribution between 8 and 12 µm.
In the present scheme, the retrieved state vector includes (T s , ) alone, whereas the principal atmospheric parameters are obtained from the ECMWF analysis and are not further iterated within the retrieval scheme.Thus, the question  arises concerning the potential bias on the retrieved parameters (T s , ) resulting from the uncertainty of those not retrieved, that is (T , Q, O).
We stress that in our scheme the (non-retrieved) atmospheric state vector (T , Q, O) is obtained from the spacetime collocated ECMWF analysis, which, especially for arid regions such as that analyzed in this paper, could be significantly in error in daytime (Masiello and Serio, 2013).
The assessment of the bias on the retrieved pair (T s , ) which arises from a non-perfect knowledge of the atmospheric state vector, can be performed through a linear perturbation analysis by dealing with a generic atmospheric parameter, say X, as a interfering factor.
Within the context of optimal estimation, (e.g., Rodgers, 2000), which (as shown in Sect.3.2.1)applies to any iteration step of the Kalman filter methodology, the sensitivity of the retrieved vector, v, to a difference X = X − X o of the given atmospheric parameter, X, with respect to the reference state X o assumed in the forward model calculations can be computed according to (Carissimo et al., 2005) www where K is the Jacobian matrix of the retrieved vector and K X is the Jacobian matrix of the interfering factor computed at the reference state X o .Equation ( 29) can be used to check the impact of possible biases in the ECMWF analysis on the retrieval for surface emissivity and temperature.To obtain realistic situations we have used a couple of day-night IASI spectra (see Fig. 10) recorded on 10 July 2010 over the Sahara desert at two close locations which are included in the target area shown in Fig. 2.These two IASI spectra have been inverted for (T s , , T , Q, O) using the so-called ϕ-IASI package (Amato et al., 1995;Masiello and Serio, 2004;Carissimo et al., 2005;Grieco et al., 2007;Masiello et al., 2009;Masiello and Serio, 2013).The IASI retrieved atmospheric state vector is compared to the ECMWF reference state vector in Fig. 10.We see that large differences arise in daytime, mostly concerning the lower troposphere.For nighttime we have a good agreement for the surface temperature (303.9K of ECMWF and 303 K of IASI), whereas for daytime we have a disagreement which is as large as 12 K (310.1 of and 321.7 of IASI).
We can take the difference, X IASI −X ECMWF , as a realistic departure of the ECMWF analysis from the true atmospheric state vector and compute, through Eq. ( 29), the resulting bias over the retrieved surface emissivity and temperature.In doing so, we have used S a defined according to with S e obtained from the UW/BFEMIS database.It should be stressed that S a defined in Eq. ( 30) gives the less favorable situation.In fact, as iterations evolve, the matrix S a evolves as well according to Eq. ( 16) and its norm tends to decrease.
In this situation, it can be shown (Carissimo et al., 2005) that the interfering effect also tends to decrease.Thus the calculations we are going to show should be viewed as an upper boundary to the impact of interfering atmospheric factors.With this in mind, Table 3 shows the impact over the retrieval of the interfering factors, (T , Q, O).It can be seen that even in this least favorable case, the impact is modest and much lower than the precision of the retrieval.As expected, the impact is larger during daytime, although mostly affecting the second decimal digit for skin temperature and below the fourth decimal digit for emissivity.
Based on this result, we have implemented the (T s , )version of the Kalman filter methodology by considering the simultaneous use of channels 4, 6, and 7.The retrieved state vectors are constructed from surface temperature and emissivity alone.We do consider atmospheric parameters in the state vector.These come from space-time collocated ECMWF analysis; however, they are not retrieved.
It is worth mentioning that the conclusion reached in this section applies to our retrieval scheme and does not have a general validity.The impact of possible interfering factors depends on the regularization determined by the matrix S a , and tends to attain its largest value in the limit S −1 a → 0 (Carissimo et al., 2005), that is, for the case of unconstrained least squares.

Assessing the performance of the emissivity/temperature retrieval
We begin with the description of results obtained by processing SEVIRI data at the three small test areas shown in Fig. 2. With this first series of results, we want to address issues such as the precision of the method and its convergence properties.
Results obtained from the analysis of the full target area will be shown later in this section.The Kalman filter has been applied to SEVIRI atmospheric channels alone.These are the channels at 12, 10.8 and 8.7 µm.The corresponding channel emissivity and the surface temperature have been simultaneously retrieved with a time resolution of 15 min.This time step also coincides with the sequential updating rate of the filter.
For the sake of clearness, Table 4 summarizes the many settings of the filter.Note that in computing background vector and the related covariance matrix from the UW/BFEMIS database, we have not used the data for July 2010, which are used for comparison with our results.day and one single SEVIRI pixel from the Sahara desert test area, and therefore corresponds to the highest space-time resolution of the methodology.It is possible to see that even for a time resolution of 15 min, temperature is obtained with a precision of ≈ ±0.2 K and better, whereas emissivity is obtained with a precision better than ≈ ±0.005.The emissivity retrieval shown in Fig. 11 corresponds to the channel at 8.7 µm.This channel is in the middle of the quartz reststrahlen band and has the higher contrast in the atmospheric window.From Fig. 11, we see that the emissivity tends to follow the daily cycle, with larger values obtained during nighttime/early morning.This is better evidenced from the analysis of the long sequence of clear-sky days shown in Figs. 12 and 13.The analysis refers to a single SEVIRI pixel in the Sahara desert (30.66 • N, 5.56 • E) and has been performed with a time resolution of 15 min.According to Li et al. (2012), we have that the amplitude of the daily cycle is quite evident for the SE-VIRI channel at 8.7 µm where peak-to-peak variations can reach ≈ 0.03.Smaller variations are found at 10.8 µm (below 0.01).At 12 µm the variation is much less pronounced and in some cases it seems to have a reverse sign with respect to the pattern at 8.7 µm, an effect which has been reported also by Li et al. (2012).
These daily variations of emissivity over desert sand, even in the dry (non-rainy) season, have been first reported by Li et al. (2012).It is very likely that these variations are the result of day-night sand evotranspiration, which occurs for direct adsorption of water vapor at the surface (Agam andBerliner, 2004, 2006;Mira et al., 2007).Clear-sky daily variation of emissivity is more pronounced for desert sand because of the strong contrast of quartz absorption band at 8.6 µm.However, it should be noted that the a posteriori covariance matrix of the analysis, Ŝt (see Eq. 11) shows a relatively strong anticorrelation of the retrieved T s and .This anticorrelation could introduce a systematic drift in the retrieved    Fig. 12. Retrieved surface temperature (bottom panel) for the Sahara desert.The retrieval has been obtained with the filter for ten consecutive days.In the legend, ECMWF T s sis refers to the surface temperature analysis at the canonic within a day, whereas ECMWF T s is the ECMWF surface t ture linearly extrapolated to the SEVIRI time steps.The upp in the figure also shows the quality of the reconstructed r (channel at 12 µm). 1 r.u.=1 W m −2 sr −1 (cm −1 ) −1 Fig. 11.Upper panel: emissivity retrieval (channel at 8.7 µm) for one day and one single SEVIRI pixel over the Sahara desert; bottom panel: same as top, but for temperature.Error bars are the square root of the corresponding diagonal elements of the covariance matrix, Ŝt .parameters, and therefore could potentially be a spurious cause of the diurnal variation seen in emissivity.Nevertheless, based on our present retrieval exercises with real and simulated observations, it seems that a persistence state model for emissivity is capable to attenuate cyclic artefacts in the retrieval (as an example, see the simulation provided in Figs. 5 and 6).Moreover, for non-arid lands we have observed situations in which the day-night emissivity variation, despite the anticorrelation of T s and , is in phase with the daily temperature cycle (e.g., see Sect.4.2), that is, the reverse of the situation we have observed for conditions over desert sand.Furthermore, IASI observations (Masiello et al., 2013) confirm that the daily variation of emissivity is a genuine feature in the data.Finally, Hulley et al. (2010) has shown that emissivity retrieval from satellite observations is sensitive to changes in soil moisture.
Figures 11 to 13 are meant to exemplify that our physical scheme is sensitive to these day-night emissivity variations.An in-depth assessment of this effect is ongoing.As already stressed in section 1, the present study mostly focuses on the novel aspect of the methodology and a comparison of its results with in situ data and other similar satellite products.
Bearing this in mind, we go back to Fig. 12, from which it is possible to see that a slight cloudiness affects the observations at the beginning of the second day.We do not skip these observations when performing the retrieval, therefore Fig. 12 shows that slight cloudiness does not bring the Kalman filter to an unstable state.In other words, the stability of the filter is not influenced by slight cloudiness, although this information is forward-propagated through the forecast.
However, overcast conditions that persist a long time (e.g., from t 1 to t 2 , with t 2 − t 1 15 min) can (e.g., because of rain) drastically change the emissivity at the endpoints t 1 and t 2 .Furthermore, cloudy radiances could be undetected, in which case serious overcast conditions could negatively influence the retrieval products.

Retrieval ECMWF T s ECMWF T s analysis
Fig. 12. Retrieved surface temperature (bottom panel) for a site in the Sahara desert.The retrieval has been obtained with the Kalman filter for ten consecutive days.In the legend, ECMWF T s analysis refers to the surface temperature analysis at the canonical hours within a day, whereas ECMWF T s is the ECMWF surface temperature linearly extrapolated to the SEVIRI time steps.The upper panel in the figure also shows the quality of the reconstructed radiance (channel at 12 µm). 1 r.u.=1 W m −2 sr −1 (cm −1 ) −1 Fig. 12. Retrieved surface temperature (bottom panel) for a site in the Sahara desert.The retrieval has been obtained with the Kalman filter for ten consecutive days.In the legend, ECMWF T s analysis refers to the surface temperature analysis at the canonical hours within a day, whereas ECMWF T s is the ECMWF surface temperature linearly extrapolated to the SEVIRI time steps.The upper panel in the figure also shows the quality of the reconstructed radiance (channel at 12 µm). 1 r.u.= 1 W m −2 sr −1 (cm −1 ) −1 .and, furthermore, only retrieval which satisfies the cost function condition of Eq. ( 13) are propagated through the filter.
To this end, it should be stressed that KF does not need to deal with equally spaced times.This is exemplified in Fig. 14, which shows the retrieval for surface temperature corresponding to whole month of July for the Seville test area.The analysis has been performed only for clear-sky soundings (according to the operational SEVIRI cloud mask) and has been spatially averaged over the 187 available SEVIRI pixels.Cloudy radiances are skipped in the analysis, which means that we use a time step which is not a constant.Missing values of the surface temperature in Fig. 14 correspond to cloudy radiances.However, undetected cloudy observations could also be processed, which can drift the filter to regions which do not correspond to the cost function below the prescribed χ 2 -threshold.Therefore, retrievals are only considered and propagated ahead only in case the cost function χ 2 = 2S has been reduced below the χ 2 -threshold.These good retrievals are shown in Fig. 14.
The retrieval for emissivity is shown again in Fig. 14 (bottom panel).Also in this case, the results have been averaged over spatially adjacent clear-sky pixels.We stress that clear sky is defined according to the SEVIRI cloud mask, which can still contain undetected cloudiness.These undetected clouds cause the occasional spikes seen in Fig. 14.It is also interesting to note, that we also see a cyclic emissivity behavior for Seville, although now the amplitude of these variations is confined within ±0.01.
Figure 15 exemplifies the analysis for the case of sea surface.The retrieval for (T s , ) refers to the Mediterranean  target area shown in Fig. 2. Also in this case the results have been spatially averaged.Possible gaps in the time sequence correspond to time intervals characterized by the presence of cloudiness.

Comparison with ECMWF T s , AVHRR-AMSR and in situ land surface temperature observations
Figure 12 shows that the ECMWF model compares fairly well with the retrieval at nighttime hours, whereas during daytime ECMWF surface temperature is biased significantly low.This is in line with the deficiencies in ECMWF model skin temperature identified by Trigo and Viterbo (2003).
To have a better assessment of this bias, we have spatially averaged the data over the ECMWF grid box of 0.5 × 0.5 degrees.In this way the results are much more consistent with the horizontal spatial resolution of the ECMWF analysis.
The results are shown in Fig. 16.For the desert site, we find that the bias at midday reaches about 9 • C and has a minimum at midnight, when the bias is about 1 • C. At 06:00 and 18:00 UTC the bias is still negative and has a magnitude of about 2 • C.
For the case of the test site of Seville, we have observed a negative bias of ≈ 7 • C at midday.However, for the other three canonical hours of the ECMWF analysis, the bias is below 1 • C.
A much better comparison has been obtained in the case of the ocean site, also shown in Fig. 16.In this case, the overall bias is about −0.3 • C (ECMWF is slightly warmer than KF).However, we also see a dependence on the hour of the day.The bias is almost zero at 00:00 and 18:00 UTC, and reaches ≈ −0.6 • C at 12:00 UTC.Thus, it seems that the ECMWF model also has a bias for sea surface, which depends on the daily cycle.
A very good agreement has also been found with the AMSR+AVHRR OISST analysis (see Fig. 17).The analysis has been used to compute the skin temperature over the small Mediterranean target area shown in Fig. 2. The results shows that SEVIRI KF captures the correct day-to-day variations of the skin temperature.Daily average temperatures agree within 0.5 • C and, whereas the agreement of the monthly average temperature is within 0.1 • C (23.93 • C SEVIRI KF vs. 24.06• C of AMSR+AVHRR OISST).The relatively large difference (about 0.7 • C) around days 23-24 is likely due to the effect of cloudiness combined to the coarser spatial resolution of the OISST analysis.Finally, we will now show and discuss the comparison with in situ land surface temperature observations at the Evora station (southern Portugal).For this station the SE-VIRI KF analysis for temperature and emissivity was computed for all clear observations, with clear sky defined according to the SEVIRI operational cloud mask.
Figure 18 (upper panel) compares the surface temperature for three consecutive clear-sky days in July 2010.It is seen that the SEVIRI KF analysis is slightly upward biased with respect to the in situ observations both at midday and before sunrise.This behavior is also obtained for the LSA SAF T s product and can be partially explained by the heterogenous scene, even though the methodology used to process in situ observations has been designed with the aim of minimizing     this effect.When we consider the comparison with SEVIRI LSA SAF product, the midday and nighttime bias tend to be confined well below 1 • C. The position of Meteosat-9 with respect to Evora, favors the observation of sunlit surfaces.The current compositing of ground data does not include an accurate weighing of sunlit and shadowed ground fractions, which also may lead to in situ temperatures being cooler than those actually observed by SEVIRI.This is further corroborated by the comparison between SEVIRI KF analysis for temperature and the SEVIRI LSA SAF T s product; two independent methodologies produce very close values, with negligible systematic differences and standard deviation of about 0.8 • C.
Emissivity retrieval for Evora (see Fig. 18), with the highest value obtained for 8.7 µm, clearly above that obtained for 10.8 µm, is consistent with the emissivity spectra for dry grass (Seemann et al., 2008;Baldridge et al., 2009).This is in agreement with the type of landscape observed around the station during the summer, when the understorey dries out completely.
Also for Evora an emissivity wavy pattern is visible, although its amplitude is very small (within ±0.005, as exemplified in Fig. 18 (bottom panel)).LSA SAF analysis uses an almost constant emissivity at 10.8 µm, which has a value of 0.975 (black line in Fig. 18).The SEVIRI KF analysis shows that the emissivity at this channel is ≈ 0.972 and varies in phase with the daily temperature cycle.Although, as previously stated, the amplitude of this variation is of the order of few parts per thousand.The full set of SEVIRI KF temperature retrievals for the Evora station is compared to in situ observations in the scatter plot of Fig. 19.This figure confirms the presence of a positive bias of 1.10 • C in the KF analysis.Again, we stress that this difference is within the uncertainty of the comparison between in situ and satellite observations.The bias is nearly zero when we compare SEVIRI KF to SEVIRI LSA SAF.

Monthly maps
We have used the scheme to perform (T s , ) retrieval over the full target area for land surface and the results are summarized in this section.
Figure 20 shows the monthly map for the surface temperature and compares it with the equivalent map derived from  the ECMWF analysis.The comparison allows us to appreciate the high horizontal spatial resolution (0.05 • ×0.05 • ) compared to that of ECMWF, which is ten times less resolved (0.5 • × 0.5 • ).Because of the monthly average, differences tend to be lower than those seen for hourly and daily values.However, especially for the arid regions differences up to 5 K are still visible.
Figure 21 shows the monthly map of the channel emissivity at 8.7 µm.The difference with the UW/BFEMIS database for the same month and geographic region is also shown in the same figure.Differences appears to be more marked for the desert sand, where the variability is much larger because of the strong response from quartz particles.However, the agreement is generally good and no important deviations are seen.The map of the channel emissivity at 8.7 µm shows very well the details of seas of sand in the Sahara desert.These correspond to the bluest areas in the map and are characterized by the lower value of emissivity.
For the sake of brevity, the maps corresponding to the other two window channels are not shown.The comparison of the results with monthly maps from the UW/BFEMIS database for the same date and location shows that differences for these channels are normally below 0.01 (Serio et al., 2013).
Finally, it is worth mentioning that the results shown in this section have been intercompared with those obtained by IASI (infrared atmospheric sounder interferometer), for the same target area and dates, in a recent paper by Masiello and Serio (2013).The intercomparison showed that SEVIRI and IASI products for temperature agree within 1 K, whereas emissivity retrievals are found highly consistent with differences normally of the order of ≈ 0.001.

Conclusions
In this paper we have described a Kalman filter methodology and its implementation for the retrieval of surface temperature and emissivity from SEVIRI atmospheric window infrared channels.
The methodology has been applied to a case study characterized by many surface features (vegetation, cultivation, urban areas, bare soil, desert sand and sea water, to name a few).
It which allows us to separate the radiative effects of the two parameters.
The analysis performed on the basis of a case study has revealed many important features regarding the time evolution of emissivity.For desert sand we observe day-night variations which are anticorrelated with the daily temperature cycle.Conversely, for other types of surface features, it seems that there is only a very small day-night variation which tends to be correlated with the daily temperature cycle.
It has been shown that the Kalman filter can handle unevenly spaced data acquisition times; this allows us to process long sequences of data in which cloudy observations are simply skipped.However, the effect of raining clouds can alter the emissivity and introduce sharp gradients in its time evolution that could be inconsistent with the persistence state equation and the relative large time scale assumed for this parameter.This effect could be alleviated by reinitializing the Kalman filter in presence of a big gap in the time sequence because of cloudiness.However, this is a point that has to be addressed with suitable case studies and therefore needs further investigation.
The results have been compared with several independent observations.These comparisons lead us to conclude that the scheme is accurate and can be reliably extended to the full SEVIRI disk.
The KF methodology for the retrieval of (T s , ) we have developed in this paper considers a situation in which we give (at least for surface temperature) poor confidence to the dynamical model (we assign a relatively high value for the stochastic variance).Therefore, the question of whether a dynamical model should be considered at all could be posed, we could use, e.g., Optimal Estimation (Rodgers, 2000).There are however mainly two reasons in favor of KF for the problem at hand.First, without the time constraint for emissivity (which allowed us to prescribe low error in the persistence model) we would not successfully separate emissivity from temperature (which are inherently anticorrelated because of the radiative transfer equation).By forcing emissivity to persistence we prevent its retrieved value from developing an unrealistic time dependence that is anticorrelated to surface temperature.Second, although with Optimal Estimation we could accommodate time continuity within the a priori background matrix through accumulation of the data on a given time slot; this would be at the expense of increasing the dimensionality of the retrieval system (Serio et al., 2013).In contrast, KF is much more efficient at keeping the dimensionality of the data space low.
It is also worth mentioning that, although we assign a relatively high value to the stochastic variance of the surface temperature, we yield a retrieval precision for this parameter (square root of the diagonal of the a posteriori retrieval covariance matrix) of ≈ ±0.2 • C. Once again, this result implies that the a priori information for temperature has little to no impact on the final estimate which is, therefore, largely dominated by the data.
To summarize, the inclusion through a KF of a persistence model within the retrieval scheme allows us to effectively separate emissivity from surface temperature because it attenuates the (T s -) anticorrelation statistical structure imposed by Least Squares estimation.In addition, the (fortunate) fact that the surface temperature is strongly constrained by the data allows us to relax the persistence constraint for this parameter alone and obtain an effective simultaneous retrieval.
The application we have considered in this paper focuses on time continuity and neglects possible spatial correlation of the surface temperature field.The use of a temporal constraint alone and the statistical independency of nearby T s grid points also reflects our present understanding of the emissivity and surface temperature space-time evolution.Nevertheless, we have shown that results are realistic using an uncorrelated assumption.However, it is important to stress that in our study the KF methodology has been established in its general form which applies to both spatial and temporal constraints.We think that the use of physical algorithms such as that developed in this study, once applied to the full disk, can lead to the formulation and/or improvements of statistical ensembles (such as the UW/BFEMIS database), which should allow us to get better insight into understanding spatial features of (T s , ) fields and, finally, enable us to apply a full 3-D Kalman filter.
It is fair to say that for the present (T s , ) problem we could also make use of a Kalman smoother (KS) (Rodgers, 2000).KS is a particular application of KF: a KF run forward and backward.The error analysis shown in this paper demonstrates that we can already achieve a very useful precision with a KF: for temperature, ±0.2K and for emissivity, ±0.005.In principle a KS approach could further improve the results, but the study was done with regard to an operational implementation where the additional logistics would be prohibitive.
Finally, we want to stress that although the case study developed in this study is limited to surface parameters (because of the limited information content of SEVIRI infrared channels on atmospheric parameters), the retrieval methodology has been described in its most general framework and can therefore provide guidance to its application to future instruments such as MTG-IRS.This instrument will have some 2000 spectral channels.Therefore, the data space, rather than the parameter space, will be driving the design of a product retrieval algorithm.Even with m = 2000, a 2-D Kalman filter (time × vertical) is feasible in terms of computational costs.In this respect, if we consider that the observational covariance matrix for MTG-IRS is expected to be nearly diagonal (which implies conditional independence of the observations), the Kalman filter update can be done sequentially (Nychka and Anderson, 2010;Rodgers, 2000).With this approach we need to store only m = 2000 diagonal elements Atmos.Meas.Tech., 6, 3613-3634, 2013 www.atmos-meas-tech.net/6/3613/2013/and use a numerical algorithm which does not involve any matrix inversion.From a computational point of view, the dimensionality of the problem would be driven by the analysis covariance matrix, S a , which at this point could include also suitable spatial constraints making the methodology 4-D.However, in case of atmospheric parameters, the issue of the dynamical model becomes much more important than for the (T s , ) problem.In fact, we know that infrared observations from satellites can lack the spatial vertical resolution to resolve, e.g., temperature and moisture structures in the boundary layer and lower troposphere.One possible approach could be to use the ECMWF analysis or forecast directly as the state equation.With satellite observations every 15 min and with 3 km resolution, the present philosophy of data-driven approach with a persistence model and suitable space-time stochastic error levels could be exploited to retrieve a 4-D atmospheric field.

Fig. 2 .
Fig. 2. Target area (rectangled in blue) used to check the retrieval algorithms.The figure also show the position of smaller test areas situated East of Seville (Spain), in a flat dune area in the Sahara desert, and below Sardinia island in the Mediterranean sea.

Fig. 2 .
Fig. 2. Target area (blue rectangle) used to check the retrieval algorithms.The figure also show the position of smaller test areas situated east of Seville (Spain), in a flat dune area in the Sahara desert, and below island of Sardinia in the Mediterranean sea.

Fig. 3 .Fig. 4 .
Fig.3.Example of overlapping between the SEVIRI fine mesh and that coarse corresponding to the ECMWF analysis.

Fig. 5 .
Fig.5.Retrieval exercise using simulations with a stochastic term for T s equal to 1 K 2 and f = 10.U temperature retrieval; lower panel: emissivity retri SEVIRI atmospheric window channels.

Fig. 3 .
Fig.3.Example of overlapping between the SEVIRI fine mesh and that coarse corresponding to the ECMWF analysis.

G
.Masiello et al.:  Kalman filter surface temperature and emissivity retrieval from geostationary platforms −5.6 −5.4 −5between the SEVIRI fine mesh and corresponding to the ECMWF analysis.
of SEVIRI atmospheric window channels 7, 6 ctively (12 µm, 10.8 µm and 8.7 µm) response (panel b) ng of the daily temperature cycle shown in panel a).Panel he derivative ratio D (see text for details) corresponding e SEVIRI atmospheric window channels.In panel b) r.u.

Fig. 5 .
Fig.5.Retrieval exercise using simulations with a variance of the stochastic term for T s equal to 1 K 2 and f = 10.Upper panel: skin temperature retrieval; lower panel: emissivity retrieval at the three SEVIRI atmospheric window channels.

Fig. 5 .
Fig. 5. Retrieval exercise using simulations with a variance of the stochastic term for T s equal to 1 K 2 and f = 10.Upper panel: skin temperature retrieval; lower panel: emissivity retrieval at the three SEVIRI atmospheric window channels.

Fig. 6 .
Fig.6.Error analysis for the retrieval exercise shown in Fig.5.The precision of the retrieval (square root of the diagonal of the covariance matrix, Ŝt ) is shown by the ±1σ tolerance interval.

Fig. 6 .
Fig. 6.Error analysis for the retrieval exercise shown in Fig. 5.The precision of the retrieval (square root of the diagonal of the covariance matrix, Ŝt ) is shown by the ±1σ tolerance interval.

Fig. 7 .
Fig. 7. Retrieval exercise similar to that shown in Fig.5, but now the variance of stochastic term for T s is equal to 0 K 2 and f = 10.

Fig. 7 .
Fig. 7.Retrieval exercise similar to that shown in Fig.5, but now the variance of stochastic term for T s is equal to 0 K 2 and f = 10.

Fig. 8 .Fig. 9 .
Fig. 8. Error analysis for the retrieval exercise similar to that shown in Fig.5, but now with the propagation operator M = 0.The variance of the stochastic term for T s is equal to 1 K 2 and f = 10.

Fig. 8 .
Fig.8.Error analysis for the retrieval exercise similar to that shown in Fig.5, but now with the propagation operator M = 0.The variance of the stochastic term for T s is equal to 1 K 2 and f = 10.

Fig. 9 .
Fig. 9. Kalman filter retrieval analysis for skin temperature as a of the stochastic variance term for T s .The retrieval has been spatially averaged over the grid box of size 0.5 × 0.5 degrees shown in Fig. 1.

Fig. 9 .
Fig. 9. Kalman filter retrieval analysis skin temperature as a function of the stochastic variance term for T s .The retrieval has been spatially averaged over the grid box of size 0.5 × 0.5 degrees shown in Fig. 1.
Fig.10.a) Day-Night pair of IASI observations over a Sahara desert site (1 r.u.=1 W m −2 sr −1 (cm −1 ) −1 ); b) temperature retrieval and comparison with the space-time collocated ECMWF analysis; c) H2O retrieval and comparison with the space-time collocated ECMWF analysis; d) O3 retrieval and comparison with the space-time collocated ECMWF analysis.

Fig. 10 .
Fig. 10.(a) Day-night pair of IASI observations over a Sahara desert site (1 r.u.= 1 W m −2 sr −1 (cm −1 ) −1 ); (b) temperature retrieval and comparison with the space-time collocated ECMWF analysis; (c) H 2 O retrieval and comparison with the space-time collocated ECMWF analysis; (d) O 3 retrieval and comparison with the space-time collocated ECMWF analysis.

Fig. 11 .
Fig. 11.Top: emissivity retrieval (channel at 8.7 µm) for one day and one single SEVIRI pixel over the Sahara desert; bottom panel: same as top, but for temperature.Error bars are the square root of the corresponding diagonal elements of the covariance matrix, Ŝt .
of the day, zero at midnight of July the first) of the day, zero at midnight of July the first) Fig. 13.(Ts, ) time evolution for the retrieval exercise shown in Fig. 12. a) Surface temperature, b) emissivity at 12 µm, c) e 10.8 µm, d) emissivity at 8.7 µm.Black dots mark the 12:00 UTC to identify the times of the emissivity minima as compared to

Fig. 14 .Fig. 15 .Fig. 14 .
Fig. 14.Retrieved surface temperature (a) and emissivity (b) for the Seville test site.Results have been averaged over the 187 adjacent SEVIRI pixels.The retrievals included are only those which correspond to clear sky soundings and χ 2 = 2S cost function values below threshold.Data have been processed (and are shown) at the SEVIRI repeat time of 15 min.

Fig. 14 .Fig. 15 .
Fig. 14.Retrieved surface temperature (a) and emissivity (b) for the Seville test site.Results have been averaged over the 187 adjacent SEVIRI pixels.The retrievals included are only those which correspond to clear sky soundings and χ 2 = 2S cost function values below threshold.Data have been processed (and are shown) at the SEVIRI repeat time of 15 min.

Fig. 15 .
Fig. 15.Retrieved surface temperature (a) and emissivity (b) for the Mediterranean Sea.Results have been averaged over the 178 adjacent SEVIRI pixels.The retrievals included are only those which correspond to clear-sky soundings and χ 2 = 2S cost function values below threshold.Data have been processed (and are shown) at the SEVIRI repeat time of 15 min.

Fig. 16 .
Fig. 16.Example of scatter plots of retrieved and ECMWF T s for the three test areas.To be properly compared with ECMWF products, retrievals have been spatially averaged over the 0.5 × 0.5 grid boxes shown in Fig. 2. a) Sahara desert; b) Seville site; c) Mediterranean basin.
Fig. 17.Comparison of the daily average sea surface temperature retrieved with SEVIRI and that computed on the basis of AMSR+AVHRR OISST analysis.Each tiny red line corresponds to one single SEVIRI pixel (178 pixels in total).The two tick lines correspond to results that have been spatially averaged over the Mediterranean target area shown in Fig.2

Fig. 16 .
Fig. 16.Example of scatter plots of retrieved and ECMWF T s for the three test areas.To be properly compared with ECMWF products, retrievals have been spatially averaged over the 0.5 × 0.5 grid boxes shown in Fig. 2. (a) Sahara desert; (b) Seville site; (c) Mediterranean Basin.

Fig. 16 .
Fig. 16.Example of scatter plots of retrieved and ECMWF T s for the three test areas.To be properly compared with ECMWF products, retrievals have been spatially averaged over the 0.5 × 0.5 grid boxes shown in Fig.2.a) Sahara desert; b) Seville site; c) Mediterranean basin.

Fig. 17 .
Fig. 17.Comparison of the daily average sea-surface temperature retrieved with SEVIRI KF and that computed on the basis of AMSR+AVHRR OISST analysis.Each tiny red line corresponds to one single SEVIRI pixel (178 pixels in total).The two tick lines correspond to results that have been spatially averaged over the Mediterranean target area shown in Fig.2.

Fig. 18 .Fig. 19 .
Fig. 18.Upper panel, land surface temperature for three consecutive clear sky days in July 2010 estimated according to this work (SEVIRI KF), in situ observations, and SEVIRI LSA SAF analysis.Middle panel, difference between retrieval and in situ observations (retrieval-in situ).Bottom panel, KF emissivity at the three SE-VIRI atmospheric window channels.For comparison the figure also shows the emissivity at 10.8 µm assumed by the SEVIRI LSA SAF analysis.

Fig. 18 .
Fig. 18.Upper panel: land surface temperature for three consecutive clear-sky days in July 2010 estimated according to this work (SEVIRI KF), in situ observations, and SEVIRI LSA SAF analysis.Middle panel, difference between retrieval and in situ observations (Retrieved-In Situ).Bottom panel: KF emissivity at the three SE-VIRI atmospheric window channels.For comparison the figure also shows the emissivity at 10.8 µm assumed by the SEVIRI LSA SAF analysis.

Fig. 19 .
Fig. 19.Upper panel: scatter plot of Ts estimated according to this work and in situ observations at the Evora station.Bottom panel: scatter plot of T s estimated according to this work and the SEVIRI LSA SAF

Fig. 20 .
Fig. 20.Left, SEVIRI monthly map of surface temperature.Right, ECMWF monthly map of surface temperature.
has been shown that by properly tuning the parameters of the state equation, we can model the different time scales of emissivity and temperature and hence develop a method www.atmos-meas-tech.net/6

Table 1 .
Definition of SEVIRI infrared channels and radiometric noise in noise equivalent difference temperature (NEDT) at a scene temperature of 280 K.

Table 2 .
Example of the matrix S e for a SEVIRI pixel corresponding to a desert site (30.66 • N, 5.56 • E).The covariance matrix has been computed for the SEVIRI channels 2 to 8 in Table1and makes reference to the emissivity vector ordered from longest to shortest wave.The element (5,5) corresponds to the channel at 8.7 µm, which is in the middle of the quartz reststrahlen band and hence is characterized by the strongest variability.

Table 3 .
Potential bias affecting the retrieval of surface emissivity and temperature due to atmospheric parameters.The bias is dimensionless for emissivity and in K for the surface temperature.

Table 4 .
Summary of the settings for the KF scheme.Surface temp.stochastic variance S ηT s As S T s for land, S T s /f 2 , with f = 10 for sea surface Observational covariance matrix, S ε Masiello et al.: Kalman filter surface temperature and emissivity retrieval from geostationary pla To avoid this effect, cloudy radiances (if detected) are just skipped within the KF scheme Atmos.Meas.Tech., 6, 3613-3634, 2013 www.atmos-meas-tech.net/6/3613/2013/

www.atmos-meas-tech.net/6/3613/2013/ Atmos
. Meas.Tech., 6, 3613-3634, 2013 Upper panel, scatter plot of Ts estimated according this work and in situ observations at Evora station.Bottom panel, scatter plot of Ts estimated according to this work and the SEVIRI LSA LST analysis.