Estimation of Land Surface Temperature under Cloudy Skies Using Combined Diurnal Solar Radiation and Surface Temperature Evolution

Land surface temperature (LST) is a key parameter in the interaction of the land-atmosphere system. However, clouds affect the retrieval of LST data from thermal-infrared remote sensing data. Thus, it is important to determine a method for estimating LSTs at times when the sky is overcast. Based on a one-dimensional heat transfer equation and on the evolution of daily temperatures and net shortwave solar radiation (NSSR), a new method for estimating LSTs under cloudy skies (Tcloud) from diurnal NSSR and surface temperatures is proposed. Validation is performed against in situ measurements that were obtained at the ChangWu ecosystem experimental station in China. The results show that the root-mean-square error (RMSE) between the actual and estimated LSTs is as large as 1.23 K for cloudy data. A sensitivity analysis to the errors in the estimated LST under clear skies (Tclear) and in the estimated NSSR reveals that the RMSE of the obtained Tcloud is less than 1.5 K after adding a 0.5 K bias to the actual Tclear and 10 percent NSSR errors to the actual NSSR. Tcloud is estimated by the proposed method using Tclear and NSSR products of MSG-SEVIRI for southern Europe. The results indicate that the new algorithm is practical for retrieving the LST under cloudy sky conditions, although some uncertainty exists. Notably, the approach can only be used during the daytime due to the assumption of the variation in LST caused by variations in insolation. Further, if there are less than six Tclear observations on any given day, the method cannot be used. OPEN ACCESS Remote Sens. 2015, 7 906


Introduction
Land surface temperature (LST) is an important variable in many research areas, such as global climate change, retrieval of soil moisture content, and ground flux [1,2].Remote sensing is an effective way of providing LSTs at the regional and global scales.Remote sensors with infrared channels, such as the Moderate-resolution imaging Spectrometer (MODIS) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), play important roles in estimating LST.Currently, four major TIR-based LST retrieval methods have been developed, including the single-channel method, the split-window method, the multi-angle method and the multi-temporal method; the accuracy of most methods is within 1 K [3,4].
Most studies focus on cloud-free conditions.However, on a regional scale, the actual weather is cloudy in most regions for more than half of the year [5].The presence of clouds significantly modifies the surface energy budget [6].Many land surface temperature products from existing satellite remote sensing data only detect cloud pixels without recording the pixel's temperature; thus, LST data are lacking in cloudy areas.However, many studies concerning the application of these data (such as drought monitoring, vegetation growth and crop yield estimation) need the entire spatial distribution of the LST over a region.Therefore, LST estimates under cloudy skies (Tcloud) are urgently needed.Because of the absorption of surface emission by clouds, Tcloud cannot be calculated directly from remotely sensed thermal-infrared information.Presently, only microwave remote sensing can be used to obtain Tcloud because it is able to penetrate clouds.Jia et al. [7] retrieved LST data based on passive-microwave remotely sensed data and achieved an accuracy of 3 K relative to the MODIS LST product.However, these methods of using microwave observations are limited because microwave remote sensing is very sensitive to surface roughness and surface moisture.
On the basis of the surface energy balance, Jin et al. [6] proposed a "neighboring-pixel" approach to estimate the LST of cloudy pixels from polar-orbiting satellite data, in which the LST of a cloudy pixel is interpolated from LST observations of surrounding clear sky (Tclear) pixels within 100-300 Km or within two days.The method is limited if the clear and cloudy pixels are not homogeneous and the atmospheric conditions are non-uniform.To overcome this weakness, Lu et al. [8] calculated Tcloud, which is interpolated from temporal-based neighboring-pixel Tclear observations and is compared using the spatial-based neighboring-pixel method.The result shows that the temporal "neighboring-pixel" method is better than a spatial approach, and the absolute error is within 1.5 K.However, one disadvantage of this approach is inevitable.Specifically, Tcloud is interpolated from Tclear of temporally neighboring pixels, while the difference in net solar shortwave radiation (NSSR) in the proposed method is obtained from spatially neighboring pixels.
The objective of this study is to estimate the LST under cloudy skies from multi-temporal remote sensing observations.The methodology is presented in Section 2, and the data, including satellite and ground-based measurements, are described in Section 3. The results and discussions of the proposed method are presented in Section 4. Section 5 provides an example of estimating LST under overcast conditions from Meteosat Second Generation (MSG)/SEVIRI data.Finally, the conclusion is presented in Section 6.

Method
Assuming the 1D periodic heating of a uniform half-space of constant thermal properties, the temperature obeys the diffusion equation: where k is thermal diffusivity (m 2 •s −1 ) and T(z,t) represents the soil temperature at a distance (z) below the surface at time (t).Under a specific set of initial and boundary conditions, the model gives the surface-subsurface temperature profiles as a function of time.Initial conditions are specified 24 h or more before the time to weaken the dependence of the model results on initial values.The lower boundary condition is that the temperature is constant at a depth of 50 cm and the upper boundary condition is the energy balance equation which is listed as Equation ( 2) [9].
where NSSR is the net shortwave solar radiation (W•m −2 ), Lnet is the net longwave radiation from ground and air (W•m −2 ), H is sensible heat transfer between ground and air (W•m −2 ), LE is the latent heat transfer between ground and air (W•m −2 ) and G is heat conducted in the soil or rock unit (W•m −2 ).A solution to Equation (1) using the cosine function is 0 ( , ) cos( ( where δ = 2 / is the damping depth of the diurnal temperature wave, td is the time at which the surface temperature (z = 0) reaches its maximum and w is the angular diurnal frequency of surface temperature in a period (which is nearly π/DD (DD is the duration of daytime)).
Combining Equation (3), the ground G (z = 0) heat transport in Equation (2), which is defined as positive in the downward direction, can be written as where ) and the other parameters are the same as those in Equation (3).
Equation ( 2) can be also written as: where S is the solar constant, τ is the transmission, A is the broadband albedo, Zn is the solar zenith, Rn is the net radiation (W•m −2 ), εs is the surface emissivity, Latm↓ is the atmospheric downward radiation, λ is the latitude, δs is the solar declination angle, NSSR is the net shortwave solar radiation (W•m −2 ), w1 is the angular diurnal frequency of solar radiation in a period that is nearly π/DD (DD is the duration of daytime), ts is the time (local time 12 h in general) at which the NSSR is maximized (NSSRmax).According to Equation ( 5), NSSR can be expressed as:  6), if A and τ are stable over a day, then a harmonic variation (cosine function) can be used to describe the diurnal NSSR evolution exactly.Note, the surface is not a Lambertian reflector, i.e., the albedo is related to the solar zenith angle, and the transmittance (τ) is not constant over a day.However, A and τ are relative stable to incoming solar radiation, and diurnal NSSR in the daytime can be only approximately described with four parameters (Smin, Smax, w1 and ts), we defined the model as diurnal solar cycle model (DSC).
According to Equation ( 3), the surface temperature can be written as Many researchers think the cosine function can only be used to describe the daytime LST, and an exponential attenuation function can be employed to fit the nighttime LST [10][11][12]; thus, the diurnal evolution of the LST (DTC) can be written as where Moreover, and T0 are the two parameters to be defined: approximately represents the minimum temperature (Tmin), while T0 represents the amplitude (approximately Tmax − Tmin, where Tmax is the maximum daytime LST).Furthermore, trs is the starting time of the attenuation (near sunset), β is the decay coefficient during the nighttime, td and w are the same as those in Equation (3).
Assuming that Latm↓ − H -LE = aT + b and that the surface longwave radiation function is linearized in the vicinity (Ti), the following function can be derived by simple mathematical manipulation of Equations ( 4)-( 6): Assuming w1 is the same as w in Equation ( 9), i.e., w1 = w, and substituting Equation ( 7) into Equation ( 9), the daytime surface temperature (T) can be expressed as follows: The method used to estimate the thermal inertia can be obtained by combining tg(w × (td − ts)) with T0 as follows: If the time series of the surface temperature and NSSR can be derived (more than six observations), three parameters (Smax, ts and w) in Equation ( 6) and two parameters (T0 and td) in Equation ( 8) can be estimated by fitting Equations ( 6) and ( 8) using least square method, then P can be obtained using w × (td − ts), T0 and Smax with Equation (12).
Thermal inertia represents the resistance to a temperature change in the upper few centimeters of the surface throughout the day, and it is independent of the local time, latitude, and season.A higher thermal inertia corresponds to a slower change in the temperature for the same incoming energy [13].Solar radiation is the primary type of energy that reaches the Earth's surface; thus, the surface temperature decreases when the sky is overcast.Because of the character of thermal inertia, the change in temperature lags behind the reduced incoming solar radiation.For example, the time at which the temperature is maximized (td) often lags behind the time at which the NSSR is maximized (ts).Meanwhile, the amplitude of the variation in surface temperature is also affected by thermal inertia.When the sky is cloudy, the reduced incoming solar radiation causes a drop in the LST; but, the time in temperature change will be postponed, perhaps lasts td − ts, and the amplitude of the temperature change increases incrementally (from 0 to 1).Assuming that the LST variations are caused by variations in insolation (ΔS), which is related to cloudiness, at the same time, the amplitudes of the LST variations are related to thermal inertia (P).Considering ΔS is less than 500 W•m −2 , and P changes from 400 to 4000 W•s 1/2 •m −2 •k −1 in most situations, Lu et al. [8] also pointed out the ratio of solar radiation change and temperature change is between 30 and 300.So, the amplitudes of the LSTs variations are enlarged by a factor of 10.The method to calculate the daytime LST under cloudy skies combining the diurnal solar radiation and surface temperature is proposed as follows: Tcloud = Tclear -10 × ΔS/P (13) where Tcloud is the LST under cloudy skies, Tclear is the LST under cloud-free skies (estimated using Equation ( 8)), and ΔS is the difference in the NSSR values between clear or cloudy skies.Considering the lag between LST variations and insolation variations, ΔS cannot attain a maximum once a cloud appears.Instead, the value reaches a maximum at some time after a cloud appears (td − ts).Moreover, the degree of the influence increases incrementally (from 0-1).Therefore, ΔS is given as follows: where Sactual is the actual NSSR at t time, Sfit is fitted using Equation ( 6) which means NSSR under cloud-free skies (NSSRclear) at t time, and t is time.
To sum up, if the time series of the surface temperature and NSSR can be derived (more than six observations), three parameters (Smax, ts and w) in Equation ( 6) and two parameters (T0 and td) in Equation ( 8) can be estimated by fitting Equations ( 6) and ( 8) using least square method, then P can be obtained using w × (td − ts), T0 and Smax with Equation (12).Meanwhile, Tclear at any time during the daytime can be estimated using Equation (8) with six parameters ( , T0, w, td, ts, and β which are inversed using more than six observations with least square method).Similarly, NSSRclear can be also obtained using Equation ( 6), and ΔS can be calculated after NSSRclear minus observed NSSR, in the end, substituting these variables (Tclear, ΔS and P) to Equation ( 13), Tcloud can be calculated.
Notably, the method of estimating thermal inertia is based on homogeneous bare soil, which is not realistic.So, the P in Equations ( 12) and ( 13) is defined as the apparent thermal inertia.Furthermore, in the process of inducing thermal inertia, the angular diurnal frequency (w) is the same for Equations ( 6) and (8).Considering the varying atmospheric conditions and underlying surfaces, w exhibits a slight difference.Therefore, the mean w (for NSSR and LST) is used to obtain the LST under cloudy skies.Figure 1 presents a flow chart for estimating LST under cloudy skies.

Data from Field Experiments
To validate the proposed method, in situ data were collected from the Chang Wu experimental station located in Shaanxi Province (107°40′E and 35°12′N), which joined the Chinese Ecosystem Research Network in 1991.The site comprises a field dominated by wheat and corn in the proximity of the meteorological station.Incoming and reflected solar radiation were measured by a pyranometer.Non-contact infrared thermometers were used to measure the surface temperature.Data for an entire year (i.e., 2012) collected at 1-h intervals were used for validation in this study.

Satellite Data
The 12 spectral channels of MSG/SEVIRI cover a range of visible, near-infrared and thermal-infrared bands.The LST under cloud-free skies can be obtained from the atmospheric windows of 3.9, 8.7, 10.8 and 12.0 µm [14][15][16].The NSSR can be estimated from the visible and near-infrared channels (0.6, 0.8, 1.6, and 3.9 µm) [17,18].In this study, the LST, Down-welling Surface Short-wave Radiation Flux (DSSF) and surface albedo products for 1 April 2012, were downloaded from the Land Surface Analysis Satellite Applications Facility [19] to determine the LST under cloudy skies.

LST under a Cloudy Sky
Field measurements collected from the Chang Wu experimental station are used to validate the proposed method; the results are discussed below.

Determining Parameters in the DTC Model
As discussed by Duan et al. [20][21][22], the DTC model can only be used to describe the LST diurnal cycle under cloud-free skies, i.e., the LST estimated using Equation ( 8) denotes the LST under cloud-free conditions.If the sky is overcast, the LST daily evolution cannot be determined using the ΔS is calculated using Equation (14) The following prior information is assumed: (1) the fitted LST is equal to or slightly higher than the actual LST; (2) the observations collected before a cloud first appears or at some time (assumed to be 2 h in this study) after a cloud disappears are given higher weights (set to 2), while other observations under cloud-free conditions are given lower weights (set to 1).
LST under cloud-free skies (at least six observations distributed throughout the day) NSSR under cloud-free skies (at least four observations distributed throughout the daytime) Six parameters are fitted using Equation (8) Four parameters are fitted using Equation ( 6) P is calculated using Equation ( 12 DTC model.Considering that the LST is the same as that under cloud-free skies before a cloud appears or sometime after a cloud disappears, we believe the DTC model parameters can be estimated using observations from cloud-free conditions.Moreover, under overcast skies, the LST is lower than the LST under clear skies [6].Therefore, some prior information can be assumed: (1) the fitted LST is equal to or slightly higher than the actual LST when inverting the six parameters; (2) the observations collected before a cloud first appears or at some time (assumed to be 2 h in this study) after a cloud disappears are given higher weights (set to 2) when inverting the six parameters, while observations under cloud-free skies are given lower weights (set to 1).In general, the six parameters can be estimated using more than six Tclear observations.However, the six Tclear observations must be distributed throughout the day, i.e., they cannot be concentrated in the morning or afternoon.Figure 2 displays the fitted results for six days under various cloud conditions.The black points denote the measured LST under cloud-free skies, and the black hollow points denote the LST fitted using Equation ( 8).The root-mean-square error (RMSE) between the measured LST and the predicted LST using Equation ( 8) are within 2 K.The fit is suitable for the daytime but not for the nighttime.Although the overcast time is different among the six days, the DTC model can be used to describe the daily LST evolution, even when most of the daytime is cloudy (e.g., 8 June in Figure 2).Prior knowledge that the fitted NSSR is equal to or slightly higher than the actual NSSR is considered when inverting the four parameters in Equation ( 6) [23].Four parameters can be obtained from more than four NSSR observations using Equation ( 6).Like the LST, the four observations cannot be concentrated in the morning or afternoon; the observations must be distributed throughout the day.The diurnal NSSR cycle is shown in Figure 3, which delineates the DSC curve for six days under various cloud cover conditions.The red hollow points denote the measured NSSR in one day, and the black points denote the NSSR fitted using Equation ( 6), i.e., the clear-sky NSSR (NSSRclear).Considering that clouds suppress the amount of solar radiation reaching the surface, it is assumed that the fitted NSSR is equal to or slightly higher than the actual NSSR.Therefore, the RMSEs between the measured and predicted NSSR are large; most RMSEs exceed 170 W•m −2 .However, the fitted curve reflects the changes in the NSSR in one day and suggests that the DSC model described in Equation ( 6) can predict the diurnal NSSR evolution on clear days.Just as shown in Figures 2 and 3, measured diurnal LST curve is smoother than measured diurnal NSSR curve, it is because LST will display the slow process of change and the time in LST change lag the time in the reduced incoming solar radiation under the influence of thermal inertia.

Estimation of LST Under a Cloudy Sky
Tcloud is calculated by the proposed method after obtaining Tclear, NSSRclear and thermal inertia using Equations ( 6), ( 8) and ( 13), and the results for the six days are shown in Figure 4.The black points denote the measured LST under cloud-free skies, the red solid points denote the measured LST under cloudy skies, and the red hollow points denote the estimated LST using the proposed method.The RMSE values of the measured and estimated LST are 1.2, 1.8, 1.05, 1.2, 0.8 and 0.72 K for April 17, June 8, June 15, June 22, August 21 and September 28, respectively.Compared with the other five days, the overcast time on 8 June was much longer, which influenced the fit (Equation ( 8)) shown in Figure 2. As a result, the absolute errors of the measured LST and the predicted LST using Equation ( 13) on 8 June are higher than those on the other five days.Compared to the RMSEs between the measured and estimated LSTs under cloud-free conditions (Figure 2), the RMSE between the measured and estimated LSTs is smaller on 28 September.Because errors between the measured and estimated LSTs under cloud-free conditions are primarily produced at night (after midnight) and the RMSE between the measured and estimated daytime LSTs under cloud-free conditions is only 1.10 K, the RMSE is smaller on 28 September.The curve of the LST predicted by Equation ( 13) is very similar to that of the measured LST under cloudy skies on 17 April, 15 June, 22 June, 21 August and 28 September; therefore, Equation ( 13) can be used to describe the response of the LST changes to variations in the NSSR.Prior knowledge that the fitted LST is equal to or slightly higher than the actual LST when estimating the six parameters in Equation ( 8) suggests that a few of the measured LSTs under cloud-free skies are less than the estimated LSTs.

Sensitivity to Errors of the Estimated LST Under Cloud-Free Skies
The LST under cloudy skies (Tcloud) is retrieved using Tclear observations in the proposed method.In general, Tclear is obtained using infra-thermal band data, and errors are unavoidable due to the uncertainty in the emissivity, unknown atmospheric conditions and the inversed method, etc.Thus, errors in Tclear are fully translated into errors in the estimation of Tcloud.Regarding the errors in Tcloud from Tclear, a sensitivity analysis is performed by adding ±0.25 K and ±0.5 K to the actual Tclear.Figure 6 shows the resulting error histograms.The RMSE values of the measured and estimated Tcloud are 1.34, 1.39, 1.36 and 1.44 K with a mean error of 0.35 K, 0.12 K, −0.57K and −0.34 K, respectively when adding −0.25, −0.5, 0.25 and 0.5 K biases to the real Tclear.We find that most errors in Tcloud are within ±2.5 K from the histogram of errors.The mean error is positive when Tclear is underestimated, while the mean error is negative when Tclear is overestimated.

Sensitivity to Errors in the Estimated NSSR
In the algorithm, the daily NSSR evolution is needed.The NSSR values estimated by the method will produce some errors, and the errors will be translated into errors in the estimation of Tcloud.Regarding the errors in Tcloud introduced by NSSR, a sensitivity analysis is performed after adding ±5 and ±10 percent NSSR biases to the actual NSSR, respectively.Figure 7 shows the corresponding histograms of errors.The RMSE values of the measured and estimated Tcloud are 1.484, 1.424, 1.42 and 1.31 K with a mean error of −0.05 K, −0.05 K, −0.13 K and −0.14 K respectively after adding −5, −10, 5 and 10 percent NSSR biases to the actual NSSR.

Application to Actual MSG-SEVIRI Satellite Data
The objective of the present work is to estimate the LST under cloudy skies from MSG data.Figure 8 provides an example of the LST retrieval using geospatial coverage between 35°N-50°N and 20°E-30°E during the MSG scan on 1 April 2012 at 11:00 UTC.The LSTs under cloudy skies are not validated due to the lack of in situ measurements and remotely sensed data from other sources.Validation of LSTs at the satellite pixel scale is challenging because LSTs can vary significantly within a pixel and change within relatively short time periods [21].However, the LSTs are generally similar to neighboring pixels.Due to the appearance of clouds, LSTs under cloudy skies are lower than that under cloud-free skies.Figure 8 presents the results from 44°-45°N and 28°-29°E, which means that the method can be used to estimate LSTs under cloudy skies.Considering the absent of LST product under cloudy conditions, comparison of LST product from LSASAF and LST estimated using proposed method under cloud-free skies is performed and the RMSE is 1.3 K with a mean error of −0.5 K.

Conclusions
With more and more geostationary meteorological satellites in operation, the estimation of LST for cloudy pixels is required.In this paper, we proposed a new method to estimate LST under cloudy skies by combining diurnal solar radiation with surface temperatures.
The RMSE between the in situ LST from the Chang Wu experimental station located in Shaanxi Province (107°40′E and 35°12′N) and the predicted value is 1.23 K; thus, the proposed algorithm can be used to calculate LST under cloudy skies.
The proposed method requires the use of LST under cloud-free skies (Tclear) and diurnal NSSR.Considering the errors in the estimated LST under cloudy conditions (Tcloud) when using NSSR and Tclear, a sensitivity analysis regarding the uncertainty of Tclear and NSSR was also performed.The results show that the accuracy of the LST retrieval can be off by 0.11 K, 0.16 K, 0.17 K and 0.21 K, respectively, when adding a −0.25 K, −0.5 K, 0.25 K and 0.5 K bias to the actual Tclear.The effect of the uncertainty in the NSSR on the LST retrieval could be approximately 0.254, 0.194, 0.19 and 0.08 K, respectively, when adding −5, −10, 5 and 10 percent NSSR errors to the real NSSR.
Finally, Tcloud was calculated using Tclear and Down-welling Surface Short-wave Radiation Flux (DSSR) and albedo products of MSG2-SEVIRI from the Land Surface Analysis Satellite Applications Facility (LSASAF) using the new method.The new algorithm can be applied to LST data retrieved from a geo-stationary satellite during cloudy conditions, and it provides the ability to reconstruct diurnal LST cycles from geo-stationary satellite observations.These cycles are particularly useful in regions where ground-based meteorological observations are scarce.
Notably, the proposed method assumes that the variation in the LST is caused by variations in insolation (which is related to cloudiness) during the daytime.Therefore, the approach can only be used during the daytime.Furthermore, at least six Tclear values in one day are needed to fit the DTC model; therefore, the algorithm cannot be used to estimate the LST when less than six Tclear observations are available.These limitations will be addressed in future research.

Figure 1 .
Figure 1.Flow chart of estimated LST under cloudy skies.
estimated at any time in one day including the time under cloudy skies using Equation

Figure 2 .
Figure 2. Comparison of the measured daily temperature evolution with that predicted by the DTC model for six days under various cloud conditions.

Figure 3 .
Figure 3.Comparison of the measured daily NSSR evolution with that predicted by the DSC model for six days under various cloud cover conditions.

Figure 4 .
Figure 4. Comparison of the measured daily temperature evolution with that predicted by the proposed method for six days under various cloud conditions.

Figure 5
Figure5shows a comparison of the measured daytime LST and the estimated daytime LST using the proposed method under cloudy skies at the Chang Wu Ecosystem experimental station in 2012 (excluding the days in which less than six Tclear observations are available).The scatter plot shows that most points are distributed near the 1:1 line, and the RMSE is 1.23 K.The histogram indicates that most errors between the measured LST and the estimated LST using the proposed method under cloudy skies

Figure 5 .
Figure 5.Comparison of the measured daytime LST and the daytime LST predicted by the proposed method under cloudy skies at the Chang Wu Ecosystem experimental station in 2012 (left: scatter plot; right: histogram of errors of the measured and estimated LST).

Figure 6 .
Figure 6.Histogram of errors in the measured Tcloud and in the estimated Tcloud based on adding biases to Tclear ((A): add a −0.25 K bias; (B): add a −0.5 K bias; (C): add a 0.25 K bias; and (D): add a 0. 5 K bias).

Figure 7 .
Figure 7. Histogram of errors of the measured Tcloud and the estimated Tcloud after adding biases to the NSSR ((A): add a −5 percent NSSR bias to the actual NSSR; (B): add a −10 percent NSSR bias to the actual NSSR; (C): add a 5 percent NSSR bias to the actual NSSR; and (D): add a 10 percent NSSR bias to the actual NSSR).

Figure 8 (
right) is the LST product from the Land Surface Analysis Satellite Applications Facility (LSASAF), and the black region in the middle of the picture denotes invalid values due to the cloud cover.The black regions are filled with the LST predicted by the proposed method (Figure8, left).

Figure 8 .
Figure 8.Comparison of LST products from LSASAF, with the LST estimated by the proposed method on 1 April 2012 at 11:00 UTC (left: estimates from the proposed method; right: the LST products).