Spatial patterns in timing of the diurnal temperature cycle

This paper investigates the structural difference in timing of the diurnal temperature cycle (DTC) over land resulting from choice of measuring device or model framework. It is shown that the timing can be reliably estimated from temporally sparse observations acquired from a constellation of low Earth-orbiting satellites given record lengths of at least three months. Based on a year of data, the spatial patterns of mean DTC timing are compared between temperature estimates from microwave Ka-band, geostationary thermal infrared (TIR), and numerical weather prediction model output from the Global Modeling and Assimilation Office (GMAO). It is found that the spatial patterns can be explained by vegetation effects, sensing depth differences and more speculatively the orientation of orographic relief features. In absolute terms, the GMAO model puts the peak of the DTC on average at 12:50 local solar time, 23 min before TIR with a peak temperature at 13:13 (both averaged over Africa and Europe). Since TIR is the shallowest observation of the land surface, this small difference represents a structural error that possibly affects the model’s ability to assimilate observations that are closely tied to the DTC. The equivalent average timing for Ka-band is 13:44, which is influenced by the effect of increased sensing depth in desert areas. For non-desert areas, the Ka-band observations lag the TIR observations by only 15 min, which is in agreement with their respective theoretical sensing depth. The results of this comparison provide insights into the structural differences between temperature measurements and models, and can be used as a first step to account for these differences in a coherent way.


Introduction
In recent decades, Earth observation by satellite has progressed from experimental to routine methods for monitoring many aspects of the hydrological cycle over land.For example, cloud water and precipitation from satellite sensors are now routinely ingested into numerical weather prediction (NWP) models (Bauer et al., 2011).Soil moisture observations are at the point of entering into operational NWP assimilation schemes (Rosnay et al., 2011).Indirect observations of the evaporative fluxes are now informing drought monitoring (Anderson et al., 2011;Hain et al., 2011).However, one crucial parameter that is missing from this list is land surface temperature (LST).Even though it has been routinely measured since the first Earth observation satellites, and physically based retrieval schemes for the above parameters must account for LST in some way, it has yet to be successfully exploited as a stand-alone input to NWP models.This is striking since LST is tightly linked (even more so than soil moisture) to land-atmosphere fluxes that are a primary prediction goal for land models within NWP systems (Bosilovich et al., 2007).
The utilization of LST observations for hydrological studies is hampered by the fact that the relationship between different model and remote-sensing-based estimates of LST are poorly understood.This is because temperature is highly variable in time and space (both vertically and laterally).While, for example, soil moisture can be assumed diurnally stable, and an observation at midnight can readily be compared to another at 6 a.m.(Parinussa et al., 2012), it is obvious that no such assumption would hold for LST.The variability with depth is a problem when comparing temperature observations from different measurement techniques that have different (and uncertain) sensing depths (and/or T. R. H. Holmes et al.: Timing of the diurnal temperature cycle models which utilize different soil layers and/or soil thermal capacities).This variability with depth together with a high spatial variability poses a challenge for the in situ validation of LST, even though thermal infrared (TIR) measuring techniques can have high spatial resolutions (up to 1 km for global MODIS products for example).The resulting uncertainty regarding LST error characteristics is especially problematic when time-varying structural errors go undetected.For remotely sensed LST estimates the bias may depend on emissivity, viewing angle, atmospheric opacity or sensing depth.In NWP models structural errors in LST may be introduced through the parameterization of the heat capacity, layer depth or estimation of surface energy balance components.In previous attempts at assimilation of temperature observations into land surface models (Bosilovich et al., 2007;Reichle et al., 2010), time-varying bias was addressed through rescaling of observations to match the autocorrelation and/or diurnal LST properties of models.While this may be the preferred strategy in a data assimilation approach where the physics of the model has to be preserved, it squanders the opportunity to correct structural errors in the surface energy balance via comparison to satellite LST observations.This represents a missed opportunity because NWP centers are known to have diurnal biases in high-value predictions like precipitation (Dai and Trenberth, 2004).
In a comparison of temperature output from three different NWP models with in situ measurements, the depth difference between model and measurements could be resolved by considering the timing, or phase, (φ) of the diurnal temperature cycle (DTC) (Holmes et al., 2012).It was shown that based on the annual average difference between φ of measurements from different depths, the effect of 5 cm depth difference could accurately be corrected for.The logic behind this is that the difference in timing between two measurement or model systems represents the integrated effect of both depth difference and soil thermal properties.This timing difference is then assumed to be accompanied by an exponential change in amplitude according to heat flow principles (Van Wijk and de Vries, 1963).
Determining φ is relatively straightforward when the sampling frequency is much higher than the daily harmonic being sampled.This is true for NWP models, and also for observations from geostationary satellites.However, for a single satellite in low Earth orbit the sampling frequency at any location is much lower: 1-2 observations per day (depending on the swath width).For such satellites we need to combine the observations from multiple platforms in order to reliably estimate φ.This was shown in Holmes et al. (2013), where vertical polarized Ka-band observations from four platforms were combined before determining the Ka-band φ.That paper showed that Ka-band observations can be used to enhance NWP temperature output, but only if the temperature series are properly reconciled in terms of timing, amplitude, and minimum of the diurnal temperature cycle.In the near future, the potential sampling of the DTC by Ka-band sen-sors will be greatly enhanced by the constellation of satellites launched under the auspices of the Global Precipitation Measurement mission (Smith et al., 2007).In this paper we determine φ of NWP surface temperature estimates with an hourly output interval as provided by NASA's Global Modeling and Assimilation Office (GMAO).We compare this to the φ as determined from an inter-calibrated record of Ka-band brightness temperatures from five satellite platforms.As a third independent data source, we use thermal infrared (TIR) LST retrievals from the geostationary Meteosat-9 satellite (centered at a longitude of 0 • ), covering Europe and Africa.
When considered as group, these three sources (i.e., NWPbased, microwave-based and TIR-based) all provide independent information regarding LST and can theoretically be integrated together (via e.g. the assimilation of TIR and microwave LST observations into the NWP model) or used to improve physical retrievals methods for ET and soil moisture which require ancillary LST information.However, before these overarching goals can be accomplished, systematic differences between these LST data sets -particularly as they relate to φ -must be understood.This will not only support efforts to combine temperature from different sources, but may also help to better tailor a given temperature set to its function within physical retrieval models.For soil moisture remote sensing it may help to better adjust the temperature measurement to the sensing depth of the band that is actually used for the soil moisture retrieval.For precipitation, it may help efforts to improve the estimation of background emissivity (Stephens and Kummerow, 2007).And finally, a proper reconciliation of thermal and microwavebased temperature may improve evaporation retrievals, such as the Atmosphere-Land Exchange Inverse model (ALEXI; Anderson et al., 1997), that currently depend on suboptimal gap-filling when clouds prevent TIR observations.
In preparation for a global merging of temperature data, this paper presents a global analysis of difference in DTC timing between Ka-band temperature estimates, TIR-based temperature estimates and NWP model output.The results of this comparison provide insights into diurnal differences between temperature measurements and models, and can be used as a first step to account for them in a coherent way.

Theory
The time between solar noon and the time of the daily maximum is here referred to as the phase of the DTC (φ) and measured in hours.Because solar noon is the time of day when the sun is at its highest point in the sky for a given location, this local definition eliminates any longitudinal dependency but also the smaller effect of variations in day length throughout the year.Accordingly, all times denoted here are given in local solar time.For a surface temperature measurement, the average φ should be directly related to the incoming radiation, with a delay (damping) that is a function of the heat capacity of the soil or vegetation layer over the measurement depth of the sensor.When comparing two temperature measurements with the same spatial extent, the measurement depth will determine the level of damping of the diurnal temperature cycle (Van Wijk and de Vries, 1963) and the measurement with the earliest peak will represent the shallowest layer.TIR measurements have a sensing depth of about 50 µm, providing the shallowest practical measurement of LST.The timing of the maximum TIR temperature is typically reported between 60 and 90 min after solar noon (Choudhury et al., 1987;Betts and Ball, 1995;Fiebrich et al., 2003).Ka-band microwave emission has been shown to be a plausible alternative to TIR measurements, with much higher tolerance for clouds but a limited spatial resolution (Holmes et al., 2009).The sensing depth for Ka-band microwave emission, with a frequency of 37 GHz, is slightly deeper than TIR and varies with soil moisture.For most land surfaces it is assumed to be around 1 mm.Accordingly, the φ derived from Ka-band emission is expected to be slightly behind that observed using TIR.Only in very dry areas with no vegetation is the Ka-band sensing depth potentially much deeper, on the order of cm's (Ulaby et al., 1986), and the φ of Ka-band can be further delayed.To illustrate the difference Fig. 1 simulates the effect of sensing depth on Ka-band φ for a wet and a dry soil.The damping of the temperature harmonic with a period of a day can be described by a phase shift (dφ) proportional to the vertical distance (dz) divided by the damping depth (z D ).The damping depth is defined as the dz over which the amplitude of the harmonic is reduced by 63 %, and is an expression of the thermal properties of the soil: where a is the thermal diffusivity (m 2 s −1 ) and f is the frequency (s −1 ) of the harmonic.For a dry soil (a = 0.15 e −6 ) Eq. ( 1) yields an estimate of z D = 6.5 cm.Estimates of the temperature sensing depth (z S ) of 1 mm for a wet soil to 1 cm for a dry soil are given in Ulaby et al. (1986).The difference in z S for microwave Ka-band emission between a dry sandy soil and a wet soil is therefore about 9 mm.Dividing this by the dry soil z D of 6.5 cm equates to a 3 h 36 min shift in φ between the soil layers at the lower ends of these sensing depths (dφz).Because the shallower soil depths weigh more heavily in the measured emission, and they have larger diurnal amplitudes, they affect the timing of the DTC more strongly.Therefore, the shift in phase (dφ) of actual measured Kaband emission, originating from the entire soil profile, is estimated as dφz/4, or 54 min (see Fig. 1).Observational evidence of this dφ will be discussed in Sect. 5.

Materials
In this study we compare three independent land temperature measurements, two from satellite measurements and one based on a global NWP model.The satellite observations include a number of Ka-band sensors with global coverage and TIR measurements from a geostationary satellite.All sets are available for the full year 2009 and are described below.

Satellite Ka-band brightness temperature
Observations of vertical polarized Ka-band (37 GHz) brightness temperatures (T Ka, V B ) are available from several satellites.For 2009 we acquired observations from the Advanced Microwave Scanning Radiometer on EOS (AMSR-E), the Special Sensor Microwave and Imager (SSM/I) and the Tropical Rainfall Measurement Mission (TRMM) Microwave Imager (TMI), and Coriolus-WindSat.Detailed sensor specifications are listed in Table 1.The coverage of all polar orbiting satellites is global, whereas the equatorial orbit of TMI extends from 38 • N to 38 • S. Each sensor has a different spatial resolution, and the location of the center of the footprint and its azimuth orientation varies between consecutive overpasses.To combine these observations, they are binned to a 0.25 • regular global grid.This resolution was chosen based on the satellite with the coarsest resolution (SSM/I; see Table 1).The value for each grid cell is the mean of all observations with a footprint center within the boundaries of the cell.The inter-calibration of these five instruments makes use of the precessing nature of TRMM's equatorial orbit.This orbit was designed to sample the diurnal variation of tropical rainfall and results in regular overlap with all polar orbiting satellites, allowing for an inter-calibration of the satellites as described in Holmes et al. (2013).
Within the microwave spectrum, T Ka, V B is the most appropriate frequency to retrieve LST as it balances a reduced sensitivity to soil surface characteristics with a relatively high atmospheric transmissivity (Colwell et al., 1983).In Holmes et al. (2009) it was further shown that an assumption of a constant land surface emissivity can be used to obtain LST estimates from Ka-band with relatively high sensitivity, if not absolute accuracy.For the purpose of this paper no conversion to physical temperature is needed since only the timing of the DTC is analyzed here, not its amplitude.This means that for this paper the assumption of constant emissivity needs only to hold over the course of the day, because there is no effect of absolute bias in LST on the analysis.Still, the linear relationship between Ka-band and LST does potentially break down under frozen soil conditions or during precipitation events.This leads to the formulation of two conditions for the Ka-band data.To avoid frost conditions T Ka, V B must be above 260 K, a rough estimate of the Ka-band freezing point determined in Holmes et al. (2009).The spatial standard deviation of all Ka-band observations within a 0.25 degree grid box is used as an indicator for measurement uncertainty; above-normal σ Ka, V is attributed to active precipitation (Holmes et al., 2013).The second condition is therefore based on σ Ka, V ; a gridbox average is rejected if σ Ka, V is more than 1 K above the annual mean for that grid box.The Ka-band temperature set is referred to in the following as T Ka .An example of the resulting sampling of the combined T Ka is given in Fig. 2 for three days in September 2009.In the same graph the NWP and infrared resource are shown; they are described below.

NWP surface temperature
The modeled temperature data set was acquired from NASA's GMAO and their Modern Era Retrospectiveanalysis for Research and Applications (MERRA) (http:// gmao.gsfc.nasa.gov/research/merra,Rienecker et al., 2011).MERRA products are generated using Version 5.2.0 of the GEOS-5 DAS (Goddard Earth Observing System (GEOS) Data Assimilation System (DAS)) with the analysis and model output both at a spatial resolution of 0.5 • latitude by Surface processes in MERRA are based on the NASA Catchment land surface model (Ducharne et al., 2000;Koster et al., 2000).Each MERRA grid cell contains several irregularly shaped tiles, and each tile is further divided into subtiles based on their modeled hydrological state: saturated, unsaturated, and wilting.The surface temperature of a grid cell is obtained by area-weighted averaging of the surface temperatures of all sub-tiles within the grid cell.The sub-tile surface temperatures are prognostic variables of the model and represent a bulk surface layer with a small but finite heat capacity.For all vegetation classes except broadleaf evergreen trees, this bulk surface layer represents the vegetation canopy and a surface layer at the top of the soil column (effective layer depth < 1 mm).
In this study we analyzed the gridded surface temperatures over land, or T NWP .The data was regridded onto a 0.25 • regular grid by means of bilinear interpolation, and the hourly output was temporally interpolated to a 15 min temporal resolution.Figure 2 gives an example of T NWP at the original hourly resolution.

Geostationary thermal-infrared-based LST
TIR-based LST (T IR ) is an operational product of the Land Surface Analysis-Satellite Applications Facility (LSA-SAF, see http://landsaf.meteo.pt).It is generated using the splitwindow channels (10.8 and 12.0 µm) of the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) on board the geostationary Meteosat second-generation (MSG) satellite; therefore a high temporal resolution (15 min) of the data is possible (Kabsch et al., 2008;Trigo et al., 2011).
These LST data are provided on a 3 km equal-area grid.When regridding onto a 0.25 • regular grid, this results in a large over sampling of the grid box.If more than twothirds of the 3 km observations are masked out for a particular location and time, then that sampling average is rejected.This threshold increases the amount of data discarded by the cloud filter for T IR .We further limit the coverage of the METEOSAT-9 to the domain covered with an Earth incidence angle up to 78 • , avoiding artifacts at large view angles.The resulting METEOSAT-9 domain covers Africa, Europe and the Middle East.In Fig. 2 the high sampling rate of T IR is apparent, as are gaps (attributed to clouds) on day 1 and 2 of this example series.

Methods
In Holmes et al. (2012) it was shown that a relative estimate of φ can be determined by fitting a simple harmonic model to a temperature time series.This worked well enough when only φ differences between temperature sets are needed.However, the values itself are hard to interpret since the actual shape of the DTC rarely resembles a perfect harmonic shape, resulting in differences between the time of maximum temperature and the peak of the harmonic.
In this paper we adapt a more sophisticated model of the DTC as described by Göttsche and Olesen (2001).This harmonic-exponential DTC model improves the fit along the cooling limb of the DTC and allows for better leveraging of observations at all times.In addition, it yields an estimate of φ that can be more readily interpreted as the time lag between solar noon and peak temperature.The DTC model by Göttsche and Olesen (2001) was originally intended as a five-parameter model, but a subsequent addition to parameterize total (atmospheric) optical depth increased this to six (Göttsche and Olesen, 2009).A recent comparison paper (Duan et al., 2012) discussed several variants of the DTC model by Göttsche and Olesen (2001) and showed that improvements were caused by adding an additional free parameter related to the day length.In this paper we have to reduce the number of free parameters to limit the computational demands and speed conversion to a solution.Therefore we simplify the five-parameter model (Göttsche and Olesen, 2001) to have only three free parameters.A comparison study (not shown) demonstrated that this does not reduce the accuracy of the mean φ as determined over longer time series.The method is detailed below.
In Göttsche and Olesen (2001) the DTC model (T par ) is parameterized as a function of the time of maximum (t m ), day length (ω), diurnal amplitude (T a ), diurnal minimum (T 0 ), change in minimum from day to day (δT ), and start of the attenuation function (t s ): Figure 3 shows an example of T par and illustrates the definitions of its parameters.The attenuation constant k is calculated by making the first derivatives of the day-and nighttime equation equal at time t s : . (3)

T. R. H. Holmes et al.: Timing of the diurnal temperature cycle
To limit the degrees of freedom and increase conversion to a solution, in this study t s is fixed so that half of the decrease in temperature (over the cooling-down limb) is described by the exponential equation.That way, t s can be calculated from t m and ω as For a given day and temperature set, first-guess estimates of T 0 and T a are determined based on the minimum and maximum recorded observations within the 24 h period from sunrise to sunrise.Similarly, δT is initialized based on the difference between T 0 of the current and the next day, if available.Solar noon (t n ) and ω are calculated based on latitude and day of year according to general solar position calculations (Cornwall et al., 2003).Defining φ as the offset between optimized time of maximum temperature and solar noon t m = t n + φ, there are three free parameters: φ, T 0 , and T a .Of these, only T 0 and T a vary from day to day; φ is assumed constant over the data range.With these assumptions it is possible to determine φ in an iterative optimization loop that minimizes the squared errors (E) between T par and the data series.In Fig. 2 examples of T par are shown, as fitted for the three LST resources used in this manuscript.
The above-described fitting of the DTC model is possible for days where the available samples are sufficient to constrain the estimates of T 0 and T a , together with at least two further samples to anchor φ.The optimization loop is therefor only applied to days where the number of samples (N (d)) is four or more.Furthermore, in order to assure a reasonable fit with the clear-sky DTC model, we limit the analysis to days where the temperature does not dip below freezing, and where the maximum temperature is close to the mean t m .Finally, in order to have an acceptable signal-to-noise ratio we focus on days where the estimated amplitude exceeds 5 K.In summary, only days (d) where the following conditions are met are included in the optimization loop to determine the timing for each temperature set: The exact thresholds applied in these four conditions are chosen to select the most optimal days for the analysis, without overly limiting the number of usable data days in a given time period.Still, the effect of these criteria on the number of usable days results in the need for a time period of a year to generate sufficient data days to compensate for the uncertainty in the satellite observations, particularly in the tropical and boreal zones with limited diurnal amplitude.The above considerations and conditions also point out the need for five different satellites to estimate the timing of T Ka .The AMSR-E satellite is crucial to satisfy condition (3) because its observation falls close to the maximum of the diurnal at 13:30.The WindSat and SSM/I satellites help constrain the early morning minimum and the afternoon cooling limb.The TMI observations increase the number of data days by helping to satisfy conditions (1) and ( 3), and are indispensable for the inter-calibration of the sensors.
To test if the sampling of T Ka results in a bias relative to the hourly sampling of MERRA, we looked at the change in apparent timing for T NWP , when only observations at the overpass times of the Ka-band set are used.The effect of sparse sampling on the φ is averaged by latitude and shown in Fig. 4. By using MERRA model output for this sensitivity analysis we cannot test for the effect of noise in the observed data.Moreover, the diurnal harmonic of T NWP may have a different bias relative to the DTC model than T Ka .This potential mismatch in DTC shape may be responsible for the bias as shown around the Equator for the first test with no noise and imperfect model.We repeat the test after removing the mismatch in the DTC model, resulting in a perfect model with no noise.Only if this assumption is valid do we expect to have no bias.These results indicate that the uncertainty as introduced by the limited sampling frequency of the Ka-band data is small enough to measure Ka-band timing from this set of sensors, but only if the DTC model accurately represents the measured T Ka .

Results
The above-described method to determine φ is applied to the three temperature records T Ka , T NWP , and T IR (described in Sect.3) for the data year 2009.The resulting (0.25 • ) maps of φ are displayed in Fig. 5 for Europe and Africa, the spatial domain of MSG-9.
On average the T Ka peaks at 13:44 with lower values over Europe and highest values over deserts and tropical rainforests.Later values of peak temperature in deserts can be explained by the deeper sensing depth of Ka-band emission under dry soil conditions (see Sect. 2).In fact, the areas with latest φ correspond closely to sand deserts; see for example the Arabian Peninsula in Fig. 5a where the Rub'al Khali Erg shows up causing an hour delay of the Ka-band φ.This feature was earlier noted in terms of day/night difference of SSM/I channels (Prigent et al., 1999), polarization difference of AMSR-E channels (Jiménez et al., 2010), and even more comparable to the present analysis, in terms of phase difference between microwave and TIR temperature (Norouzi et al., 2012).In this later study the third component derived from a principal component analysis (PCA) was attributed to the phase of the diurnal cycle.The close resemblance of the spatial features generally supports this conclusion.
More surprisingly in Fig. 5a is the delay in φ over the tropical forest, which cannot be attributed to sensing depth  variation.Another unexpected feature in the maps of φ(T Ka ) is the earlier phase over mountainous areas (e.g., the Rif mountains in NW Africa, the Caucasus).Both features will be discussed in relation to the φ of T IR and T NWP .
For T NWP , with a mean peak temperature at 12:50, there is much less spatial variation in φ with the notable exception of tropical rainforest; see Fig. 5b.The distinctly different results over rainforest with values around 13:30 are explained by a higher heat capacity as parameterized for areas classified as tropical forest in the MERRA model.At more northern latitudes, higher φ values are also found over the forest areas of the eastern European plain.Areas with lowest φ seem to correspond with high-elevation areas.
Figure 5c shows the φ as determined based on T IR .The T IR data peak on average at 13:13.The higher φ values in the tropical zone match well with the delay as found for T NWP , although for T IR the area with delayed φ is slightly more  Figure 5d shows a north-south transect of the φ as determined for T Ka , T NWP , and T IR -all averaged over the longitude extent of 0-30 • E. All sets have a delayed φ around the Equator, which for T NWP is explained by the higher heat capacity as parameterized for tropical rainforest.For T Ka and T IR , with little to no penetration of the canopy layer, the delay in φ might more plausibly have to do with the effect of a diurnal pattern in cloudiness which can cause a delay in the peak solar radiation.If this is true, then MERRA predicts a correct delay in timing but for the wrong reason.Superimposed on this tropical signature, T Ka clearly shows a delayed φ over (seasonally) dry areas in northern and southern Africa, explained by the deeper sensing depth.
The patterns in φ and differences between the sets seem to align with particular land surface types.Therefore, MODIS land cover maps (MCD12C1, Version 051) are used to study φ by land surface type (International Geosphere-Biosphere Programme (IGBP) classification).Figure 6 lists the average values as obtained over selected surface types between longitudes 0 to 30 • E. This indeed confirms the general relation with vegetation type -all sets exhibit delayed φ over broadleaf forest compared to the cropland and Europe selections.Only for T Ka are delays in φ recorded over the savannah and desert land cover types.
As a consequence of the general agreement in φ between T NWP and T IR , the spatial map of their difference in φ is remarkably homogenous: φ(T IR ) − φ(T NWP ) = 27 min (±15 min) for 79 % of the land area (see Fig. 7a).Areas with larger differences are found at the edge of land masses and mountain ranges, where φ(T IR ) can be up to 1.5 h behind φ(T NWP ) as shown in Fig. 7. Interestingly, these areas also show up in Fig. 7b, sometimes even indicating that T Ka peaks earlier than T IR , which is not physically realistic.Because of this, we look to T IR for the explanation of these anomalies.Insofar as these features line up with mountain ranges, they are tentatively attributed to azimuth angle effects on T IR .Since they are retrieved from a geostationary satellite (MSG-9) located at the prime meridian, variations in azimuth angle might explain the earlier φ over mountains to the northwest of the satellite and the delayed φ over ranges that are to the east of the satellite.Such azimuth angle effects are muted within T Ka as this is composed of observations with varying azimuth angles (and obviously play no role in the NWP record).Along the coast in the tropical zone large negative anomalies show up in T IR relative to T Ka and T NWP that cannot be explained by mountains.In these areas there are few days without clouds, and we therefore attribute this to a failure of the cloud mask in T IR .
The timing difference between T Ka and T IR is greatest over the driest parts of Africa, where T Ka is 57 min (±25 min) behind T IR in φ; see Fig. 7b.This average dφ agrees with the theoretical calculations of Sect. 4. This suggests that the large increase in φ Ka over dry areas can be quantitatively explained by an increase in temperature sensing depth.As expected for all but the driest soils (especially when covered with vegetation) the difference is much closer to zero over the rest of Africa and Europe.For example, over Europe the φ maps show that T Ka is only 14 min (±12 min) behind T IR .
In order to test if the discussed patterns in φ are stable throughout the year, the procedure to calculate φ was repeated for 3-month seasonal periods.The seasonal results confirmed the large-scale north-south patterns as shown in Fig. 5d, giving confidence that land cover type is the main determinant for φ, rather than seasonal varying factors like soil wetness or cloudiness.Furthermore, it suggests that φ may be considered relatively constant in time; the spatial standard deviation of the seasonal anomaly from the annual φ is 6 min for T NWP , 7 min for T IR , and 14 min for T Ka .

Global results
This section will briefly discuss the spatial patterns in timing on the global scale, as calculated for T Ka and T NWP .The geostationary TIR resource is not yet available as a consolidated global data set and is therefore not included here.The global maps of φ confirm the land cover features as found over Africa and Europe and discussed in Sect. 5.In the global φ map of T Ka (see Fig. 8,top), delays of up to an hour from the mean φ of 13:37 for desert areas are matched with similar delays in the deserts of Australia and central Asia.Similarly, the apparent delay in timing over tropical forest is repeated over the Amazon.
For T NWP (Fig. 8, middle panel), φ is generally between 12:30 and 13:00 with three main exceptions.Tropical forest show up with highly delineated delays in φ of up to 40 min over the Amazon and Indonesia, in close agreement with the above discussed Congo Basin in Africa.As discussed earlier, this is explained by the increased heat capacity of the surface layer as parameterized for tropical forest by MERRA.Because this is a static feature of MERRA, the resulting spatial pattern of φ is expected to remain stable from year to year, as long as the same land classification is used.
Another feature that shows up is earlier φ over mountainous and high-elevation areas.Although this feature was discussed earlier, it is much more pronounced over the Andes and Himalaya.Since it appears in both T Ka and T NWP , it most likely reflects an actual pattern in diurnal heat exchange.Possible physical mechanisms for this include the effect of cooler air temperatures at higher elevation on sensible heat flux, diurnal patterns in orographic cloud formation, and slope effects on incoming solar radiation as described in Senkova et al. (2007).
A new feature that appears in these global maps is anomalies of up to an hour from the mean φ as calculated for T NWP , lining up closely with the boreal climate zone above 60 • N.This surprising feature is not matched in T Ka .Since our analysis is restricted to days with minimum temperatures above freezing and a diurnal amplitude of 5 K or higher, the φ calculated in these northern regions (with weak diurnal radiative forcing) is based on a relatively limited amount of long cool days.Unfortunately, TIR observations are not available at those high latitudes to positively attribute this anomaly to T NWP .On average T Ka peaks 40 min later than T NWP (Fig. 8, lower pane).In the temperate climates and the tropics the dφ is close to this average.Bigger differences appear in the desert areas where Ka-band peaks up to 2 h after T NWP , corresponding to a deeper sensing depth.The boreal areas show smaller differences resulting from the delayed φ of T NWP as discussed above.Due to the close link with land cover and geological features it is expected that these patterns of dφ will remain relatively stable from year to year.
Each of these features deserves discussion in more detailed studies where ground measurements are available for comparison.Such regional studies would benefit from relaxing some of the parameters that were fixed in this global study, allowing for improved fit of the DTC model (T par ) over particular biomes.In support of the main features that were dis-cussed in this paper, we conclude with two graphs that summarize the overall fit of T par to the three LST estimates over four key biomes; see Figs. 9 and 10.The displayed data reflect the annual average by hour of day for a single grid cell.Observational data are corrected for the effect of sparse sampling on the average.These graphs show that even though the fit of T par is not perfect at all times of day, it does not appear to impact the relative conclusions of the study.All sets have a clear delay in φ over tropical forest that is picked up by the fitted T par in the presented examples.In contrast, the late φ of T Ka stands out from the other sets in the desert example.In boreal forest T NWP has a relative delay in φ that is not reflected in T Ka .

Conclusions
This study demonstrates that the timing of the diurnal temperature cycle can be reliably determined from temporally sparse data sets.The method is applied to a yearlong record of geostationary TIR observations, model output, and a combination of low Earth-orbiting satellites with microwave radiometers.The spatial patterns in φ for each temperature set are explainable based on consideration of land surface type and basic physics describing the penetration depth of microwave observations.An interesting observation is that over tropical forest the timing is delayed by 30 to 40 min relative to the average phase in both satellite data sets.While this delay is accurately modeled in MERRA, it may be for the wrong reason if the delay is caused by a diurnal variability of cloudiness rather than an increase in heat capacity of the sampled surface layer.On the other hand, deserts cause a delay in the timing of Ka-band diurnal temperature cycle which is not matched in the TIR nor MERRA estimates.This is explained by a deeper sensing depth for the microwave observations in dry soils and becomes especially extreme in dry sand deserts.The timing of Ka-band observations outside of the desert and semi-desert areas is (on average) 15 min after TIR.This small delay of the DTC compared to TIR agrees with a slightly deeper sensing depth for Ka-band, 50 µm vs. 1 mm.On the other hand, MERRA seems to model the average peak of the DTC about 23 min before that of TIR, which should be the shallowest observation possible of the land surface.Therefore, if the goal is to model the surface temperature, it appears the heat capacity of the surface layer is set slightly too low in MERRA.To put this timing difference in context of biases in T, there are two main considerations.First, if this timing difference is indeed a result of a low bias in heat capacity, then the associated overestimation of the diurnal amplitude would be 10 % for a 23 min timing difference (Holmes et al., 2012).Secondly, the timing difference will introduce a diurnal bias term.These two effects together result in a harmonic diurnal bias with an amplitude that depends on the diurnal temperature range (for example, 1.4 K bias for a T with 20 K diurnal range).This type of time-variant structural bias terms are much harder to account for in data assimilation approaches than constant bias terms.
This study has identified structural differences in diurnal timing between MERRA, TIR and Ka-band-based land surface temperature estimates and constitutes one of the first global analyses of the effects of vegetation and sensing depth on the timing of different temperature measurements.Even though the global maps of the timing for each set are based on the year 2009, these features are expected to be relatively stable from year to year.The presented analysis of the timing of the diurnal temperature cycle offers a means to account for time-variant bias terms between temperature records in a physically consistent way.With these maps we can now reconcile temperature records in terms of their diurnal timing, opening the way for studies that look at differences in diurnal amplitude and daily minimum, and ultimately for a global merger of temperature data sets.

Fig. 1 .
Fig. 1.Simulation of effective temperature for Ka-band emission for a dry (T eff (dry)) and a wet soil T eff (wet)).Assuming the shallow and deep layers are weighted equally in T eff (dry), the resulting effective phase difference (dφ) is 54 min, about a fourth of the difference (dφz) between the temperatures at their respective damping depths (z = zs).

Fig. 2 .
Fig. 2. Three days of diurnal cycles of T Ka , T NWP , and T IR in September 2009.All variables are explained in Sect.3. Solid lines represent the fitted DTC, as described in Sect. 4.

Fig. 3 .
Fig. 3. Main characteristics of the diurnal temperature cycle and the definition of the phase of the DTC (φ).

Fig. 4 .
Fig. 4. Effect of sparse temporal sampling (at T Ka observation times) on retrieval of φ.
2. T 0 (d) above freezing point, 3. t m (d) within 2 h of mean t m , and 4. T a > 5 K.

Fig. 6 .
Fig. 6. 2009 mean φ determined for T Ka , T IR , and T NWP within selected IGBP land cover types.Between T NWP and T IR there is a fairly constant dφ of 20 min.T Ka agrees with T IR over land surface types with little or no barren surface.

Fig. 9 .
Fig. 9. Example of overall fit T par to the three LST estimates for selected locations.The displayed data reflect the annual average deviation T = T − (T 0 + T a /2) for each hour of day.

Fig. 10 .
Fig. 10.Example of overall fit of T par to LST estimate in the boreal region.The displayed data reflect the annual average deviation T = T − (T 0 + T a /2) for each hour of day.

Table 1 .
Specifications of satellite sensors providing Ka-band microwave observations used for land surface temperature estimation.