Components of near-surface energy balance derived from satellite soundings – Part 1 : Noontime net available energy

This paper introduces a relatively simple method for recovering global fields of monthly midday (13:30 LT) near-surface net available energy (the sum of the sensible and latent heat flux or the difference between the net radiation and surface heat accumulation) using satellite visible and infrared products derived from the AIRS (Atmospheric Infrared Sounder) and MODIS (MODerate Resolution Imaging Spectroradiometer) platforms. The method focuses on first specifying net surface radiation by considering its various shortwave and longwave components. This was then used in a surface energy balance equation in conjunction with satellite day–night surface temperature difference to derive 12 h discrete time estimates of surface system heat capacity and heat accumulation, leading directly to retrieval for surface net available energy. Both net radiation and net available energy estimates were evaluated against ground truth data taken from 30 terrestrial tower sites affiliated with the FLUXNET network covering 7 different biome classes. This revealed a relatively good agreement between the satellite and tower data, with a pooled root-mean-square deviation of 98 and 72 W m for monthly 13:30 LT net radiation and net available energy, respectively, although both quantities were underestimated by approximately 25 and 10 %, respectively, relative to the tower observation. Analysis of the individual shortwave and longwave components of the net radiation revealed the downwelling shortwave radiation to be main source of this systematic underestimation.


Introduction
An important manifestation of climate change is widespread alteration of the composition of the energy balance at the Earth's surface (Trenberth et al., 2009;Wild et al., 2013).Given the importance of being able to predict the consequences of climate change, both measurement and modelling of the components of surface energy balance attract significant attention from a broad range of related scientific disciplines (Stephens et al., 2012).Two such disciplines are hydrology and meteorology, which share a common interest in resolving the balance between sensible, H , and latent, λE, heat fluxes over a broad range of spatial and temporal scales (Anderson et al., 2012).
Net available energy, , is a core variable used to predict the magnitude of H and λE given that it is defined as the sum of these two fluxes (Wright et al., 1992;Migletta et al., 2009;Anderson et al., 2012), Published by Copernicus Publications on behalf of the European Geosciences Union.
= λE + H. (1) The utility of this definition arises from being able to also specify as the difference between the net broadband radiation, R N , and the rate of heat accumulation, G, below the plain across which R N is specified, Given that R N is routinely measured using net radiometers, this affords an opportunity to specify and hence either H or λE.For example, in modelling studies λE is invariably specified as a function of using the ubiquitous equations such as those of Penman (1948) for open water or Monteith (1965) for land surfaces (Mu et al., 2011;Mallick et al., 2014a).Despite being the rate of change of heat stock in terrestrial environments, G is often interpreted as the "ground heat flux", and attempts to measure this using heat flux plates are commonplace (Mayocchi and Bristow, 1995;Sauer and Horton, 2005;Heitman et al., 2010).These measurements prove somewhat less reliable than R N due to greater spatial heterogeneity in ground heat uptake (Gao et al., 1998;Tittebrand and Berger, 2009;Verhoef et al., 2012) allied to the fact that significant heat capacity resides in other elements of the land surface (Ochsner et al., 2007).As a result, G proves problematic in surface energy balance studies and is either ignored (Foken et al., 2006;Foken, 2008) or treated somewhat superficially (Choudhury, 1987), despite being significant under a broad range of conditions (Santanello and Friedl, 2003;Ochsner et al., 2007).Large-scale estimates of G are useful in the context of regional and global evapotranspiration modelling and for verification of regional and global circulation models (Kergoat et al., 2011).The arrival of satellite retrievals for many of the components of R N has opened up opportunities to develop largescale estimates of this variable and hence λE (Batra et al., 2006;Mu et al., 2007, Anderson et al., 2012).For example, retrievals for the components of R N have been available through the International Satellite Cloud Climatology Project (ISCCP) (Pinker and Laszlo, 1992;Stephens et al., 2012), the Earth Radiation Budget Experiment (ERBE) (Priestley et al., 2011), and the Clouds and the Earth's Radiant Energy System (CERES) (Mlynczak et al., 2011;Chen et al., 2013) onboard NASA's Earth Observing System (EOS) and Tropical Rainfall Measuring Mission (TRMM) satellites (Wielicki et al., 1998).Several studies have reported the estimation of R N using a combination of MODIS (MODerate Resolution Imaging Spectroradiometer) atmospheric and land products over the USA, China and India (Cai et al., 2007;Mallick et al., 2009;Bisht andBras, 2010, 2011) or NOAA-14 (National Oceanic and Atmospheric Administration) data over the Tibetan Plateau (Ma et al., 2002).
Unfortunately, in the absence of direct observations of G at spatial scales and coverage of satellite R N , retrievals for have had to rely on the parameterisation of G using surface temperature, albedo and vegetation index information (Bastiaanssen et al., 1998;Batra et al., 2006) or by assigning some fixed proportion of R N (Choudhury, 1987;Humes et al., 1994) in satellite-based surface energy balance models (Mecikalski et al., 1999;Anderson et al., 2012).However studies by Murray and Verhoef (2007), Hsieh et al. (2009) and recently Verhoef et al. (2012) also demonstrated that G is, by definition, a highly dynamic quantity, and that the ratio G / R N can range anywhere from 0.05 to 0.50 depending on the time of day, soil moisture and thermal properties, and vegetation density.Therefore, methods that are able to provide defensible estimates of G in conjunction with R N would clearly be of great benefit to this area for determining directly from satellite data and without relying unduly on any offline calibration.In this paper we present a method for retrieving R N and based on exploring both satellite radiance data and day-night surface temperature difference.The approach is necessarily simple in order to avoid over-reliance on models in the pre-processing and to reflect the fact that the focus of this work is the production of satellite estimates of monthly midday (13:30; all times listed are in local time, LT) for use in a simple Bowen ratio (Bowen, 1926) λE specification framework as detailed in a companion paper by Mallick et al. (2014b) (we refer to this as M2 hereafter).Taking advantage of the extensive network of terrestrial eddy covariance tower sites (Baldocchi et al., 2001) which record direct measurements of R N , H and λE, we use these measurements to derive independent non-radiative estimates of in order to critically evaluate our satellite estimates of this quantity.
The method we present here for estimating R N contrasts with more sophisticated model-based approaches which attempt to accommodate the complexity of atmospheric radiative transfer explicitly (e.g.Fouquart and Bonnel, 1980;Mlawer et al., 1997;Bisht andBras, 2010, 2011;Hou et al., 2014).There are several reasons for adopting this stance.Firstly, the estimates of R N need to be in agreement with the simple dynamic energy balance used to accommodate G when specifying .Secondly, we believe it to be important that the complexity of the methods used here is commensurate with those methods used in the simple Bowen ratio approach as described in M2.Related to this, we have tried to restrict the approach to largely using only AIRS (Atmospheric InfraRed Sounder) data which provide the satellite soundings required for the Bowen ratio estimates.This single-platform approach is adopted to ensure the estimates do not suffer unduly from blending different data sources.Finally, as is the case with this method, complex radiative transfer approaches are also prone to the effects of uncertainty (Betts et al., 1993;Morcrette, 2002;Seidel et al., 2010), and therefore the parsimony implicit in the methods used here may be seen as advantageous.

Satellite data sets
In the present study, two different data sources were used for the estimation of R N and , AIRS and MODIS.The AIRS sounder is carried by NASA's Aqua satellite, which was launched into a Sun-synchronous low Earth orbit on 4 May 2002 as part of NASA's Earth Observing System.It gives near-global coverage twice daily at 01:30-13:30 LT from an altitude of 705 km.Level 3 standard monthly day-night data products of air temperature and relative humidity profiles, cloud cover fraction, surface emissivity, near-surface air temperature, and surface-skin temperature and columnar total precipitable water at 1 • × 1 • spatial resolution were obtained for 2003 from the online data archive of AIRS, distributed through NASA Mirador data holdings (http://mirador.gsfc.nasa.gov/).The monthly products are simply the arithmetic mean, weighted by counts, of the daily data of each grid box.The multi-day merged products have been used here because the IR retrievals are not cloud-proof and the multi-day product gave decent spatial cover in light of the missing cloudysky data.The data products were obtained in hierarchical data format (HDF4) along with their latitude-longitude projection.It is also important to mention that the daily 1 • data contain orbital gaps and cloud contamination.In the 8-day data the co-incident land surface temperature in both the day and night pass was missing, and the atmospheric soundings were also missing in many places.It is the monthly data set where the soundings as well as both the day-night land surface temperatures were available and the data have complete global coverage.
We have used the MODIS Aqua atmospheric product data sets (MYD08_D3) (http://modis-atmos.gsfc.nasa.gov/index.html) at 1 • × 1 • spatial resolution for extracting the solar zenith angle field.AIRS data do not contain any surface albedo field.For generating the surface albedo fields we used narrowband surface reflectances from combined MODIS Terra-Aqua 16-day data (MCD43C4) products acquired from the MODIS data archive (http://ladsweb.nascom.nasa.gov/data/search.html).The native spatial resolution of the MCD43C4 data sets is 0.05 • .Therefore, all the narrowband surface reflectances were first resized into 1 • by 1 • to make them compatible with the AIRS spatial resolution and then the broadband surface albedo was generated from the narrowband reflectances following Liang et al. (1999) (presented in next section).It is important to mention that MODIS global albedo product (MCD43C3) contains bi-hemispherical reflectance (white-sky albedo) and directional hemispherical reflectance (black-sky albedo).Bluesky albedo can be determined by weighting the white-and black-sky albedo with diffuse skylight fraction, which is a function of the aerosol optical depth and solar zenith angle.Look-up-table-based aerosol information and parameters are needed to convert the reflectances into the blue-sky albedo.
But there are established formulations (Liang et al., 1999;Liang et al., 2002) to directly convert the narrowband reflectances into the broadband visible albedo that does not depend on any atmospheric variables and look-up tables, and therefore narrowband surface reflectances are used in the present study.
One of the core objectives of the work is to explore the potential of atmospheric sounding data.AIRS is the only dedicated sounder available which can be explored to address the objectives in the paper.MODIS also has soundings, but it was not designed for this and only has low-quality air temperature soundings.Coarse spatial resolution of AIRS would introduce many difficulties when it comes to the evaluation, but the most important aspect of the two companion papers (we refer to the current one as M1) is to introduce the possibility of using atmospheric sounding data as a means of observing surface energy fluxes (the companion paper on latent and sensible heat flux, M2).We have restricted derivation to (largely) AIRS data (we used MODIS albedo because AIRS does not contain any albedo field) in order to exploit a single platform for the entire framework.We would also emphasise that the retrievals are on one time slot per day for 13:30 LT, which is a standard for the studies that use polarorbiting satellites.

Net radiation
The approach for estimating R N uses the AIRS radiation products, although we have also made use of the MODIS surface reflectance and solar zenith angle products where necessary.R N is generated by considering the following balance between net shortwave (R NS ) and longwave (R NL ) radiation at or near the Earth's surface, where α is the surface albedo, R L↓ and R L↑ are the downwelling and upwelling thermal radiative fluxes, and R S↓ is the downwelling shortwave radiative flux (all fluxes specified in W m −2 ).Our chosen reference level for R N is the near surface given that this corresponds to the flux-based tower estimates we used in the evaluation.Therefore, surface R S↓ was estimated from its top-of-atmosphere clear-sky counterpart R S0↓ and AIRS cloud cover fraction (f ) following Hildebrandt et al. (2007), where τ A is the clear-sky transmissivity of the atmosphere, which we assume is 0.75 (Cano et al., 1986;Thornton and Running, 1999;Hildebrandt et al., 2007;Gubler et al., 2012).
Although clearly a simplification, a constant clear-sky transmissivity is widely used (e.g.Massaquoi, 1988;Bindi et al., 1992;Choudhury, 2001;Hildebrandt et al., 2007;Mallick et al., 2009)  Eq. ( 2) should help accommodate the effects of variations in both the aerosol optical depth (Kaufman and Koran, 2006;Quass et al., 2010) and atmospheric water vapour (Adhikari et al., 2006).The terrestrial surface albedo was generated using the MODIS Aqua-Terra surface reflectances r i following Liang et al. (1999) et al., 1999, 2002).The albedo of the ocean varies according to the cosine of solar zenith angle (Jin et al., 2004).Given that the oceanic surface reflectances are not available in either MODIS or AIRS, a constant albedo of 0.04 was assumed for oceans as satellite radiances are nadir.Many of the longwave components of the radiative balance are very closely related to the raw IR radiances being measured by AIRS.Because these are not in the public domain, we have attempted to recover them as follows, although in future we would anticipate using the raw IR radiances more directly if possible.R NL was calculated as where σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ), T c is the columnar air temperature, and ε C and ε S are the column and surface emissivities.Among the different schemes for calculating ε C we have used the formulation proposed by Prata (1996) as this appears to be the most reliable (Niemela et al., 2001;Bisht andBras, 2010, 2011).This scheme uses AIRS total precipitable water (ξ ) (cm) information to estimate ε C as The columnar air temperature T C in Eq. ( 6) is taken as the average of the 2 m and 1000 hPa pressure level AIRS temperatures in an attempt to reflect a weighting toward the lower troposphere when specifying R L↓ .T S and ε S are taken directly from the AIRS skin temperature and surface emissivity products.

Surface heat capacity, ground heat flux and net available energy
The definition of G stems from consideration of the nonsteady-state surface energy balance, where c is the aggregate surface system heat capacity (MJ m −2 K −1 ).The AIRS sounder platform samples twice daily at 01:30 and 13:30 LT.Despite being somewhat coarse, if a discrete time is taken, backward difference approximation of Eq. ( 8) with a sample interval of t = 12 h equivalent to that of the AIRS pass gives where T s is the day-night surface temperature change, b 1 = t/c and b 2 = − (t) t/c.If we assume that the system is approximately in equilibrium over a 24 h cycle, and that ≈ 0 at 01:30 LT (for all 30 sites analysed in this study (01:30 LT) < 0.05 (13:30 LT); see also Tamai et al., 1998;Mamadou et al., 2014), then this gives the following simultaneous equations: which can be solved analytically to derive b 1 and b 2 and hence and c for each grid cell in the AIRS global array.Equation ( 10) is a coarse approximation of Eq. ( 8) and hence potentially suffers from a number of deficiencies.Firstly, diurnal symmetry in T s is only appropriate when one considers weekly or monthly average behaviour, and that there are no additional heat losses to or gains from stores beyond the domain defined by the single heat capacity c.In this study we examined the monthly average behaviour because AIRS only gives partial global coverage on the daily timescale due to both cloud effects and the non-overlapping swath width of the sensor.Interactions with additional longterm heat stores is an issue in systems such as the oceans, where there can be a persistent heat loss/gains to/from deeper water over timescales of weeks to months, although relative to the diurnal fluctuation of stored surface heat this tends to be small (Stramma et al., 1986).Secondly, can be either positive or negative at 01:30 LT, although it tends to be only a fraction (it never exceeds 5 % of afternoon ) of its 13:30 LT value due to the supply of relatively low amount of energy at night compared to the day.This may be less true for areas of land in the height of winter with cloud-covered days and over the sea where significant daytime heat accumulation could in part be re-released as night-time latent and sensible heat.Thirdly, Eqs.(9 and 10) attribute the magnitude of the daytime G to the night-time net longwave radiative balance, which is obviously rather uncertain.Fourthly, the air emissivity computation using Prata's equation was developed for the daytime and using it for the night-time emissivity may introduce errors.Finally, all the terms in Eq. ( 8) are highly dynamic and yet are treated as constant or varying linearly over the 12 h sample interval.It is difficult to predict what the consequences of this are, as it depends on the pattern of radiative forcing throughout the day, which can vary significantly in both time and space.Some illustrative examples of the theoretical assumptions of Eq. (10a, b) are depicted in Figs. 2, 3 and 4. Figure 2 shows the examples of the diurnal symmetry of T S for clear days in three different seasons where the saw-tooth pattern between noon (13:30 LT) and night (01:30 LT) T S is evident.This clearly shows how well these two T S samples capture the dynamic range of the day and hence the discretisation is representative of the daily energy balance.Figure 3a to d illustrate the diurnal evolution of during three different times of the year (spring, summer and winter) for four broad biome categories (grassland, cropland, forest and savanna), which clearly indicates ∼ = 0 at 01:3 LT0 and within less than 5 % of afternoon .Lastly, Fig. 4 highlights the two-dimensional relationship between the noontime G (13:30 LT) and night-time R NL (01:30 LT) for the above-mentioned four biomes and the correlation between the two varied between 0.32 and 0.60, having high correlation over grassland and savanna and moderate correlation over forest and cropland.Despite large differences in the footprint size between G and R NL measurements, the inverse relationship between the two variables in Fig. 4 clearly indicates the dependence of noontime G on the night-time longwave radiation balance.Therefore, although the theoretical approximations in Eqs. ( 8) and (10a, b) seem to be somewhat coarse and might have a few limitations (as described earlier), but Figs. 2, 3 and 4 indicate strong connectivity between T S , 01:30 LT R N (or R NL ) and 13:30 LT .However, given the structure of the atmosphere and the very small energy fluxes involved, high-latitude estimates from this method are likely to be problematic in any case.

Sensitivity analysis
A general sensitivity analysis was carried out in order to assess the effects of the propagation of uncertainty onto the estimates of G, R N and .For this analysis the input terms were assigned uniform prior distributions of ±10 % for all parameters other than temperatures for which ±1 K uniform prior distributions were assumed.These assumed ranges re- semble the stated uncertainties as given in the AIRS support literature (Aumann et al., 2003;Hearty et al., 2014).The sensitivity of each output to each input was calculated assuming an average, locally linear sensitivity.These were expressed as the change in output per unit change in input, normalised by the median value of each.Only absolute sensitivities > 0.1 were considered significant.The standard deviations of the estimated distributions of G, R N and were used as the summary statistic for the measurement uncertainty of the proposed methodology.

Evaluation of R N and
To evaluate the satellite values of R N and , we have made use of the extensive FLUXNET terrestrial tower network (Baldocchi et al., 2001).Clearly, there is a scale conflict here, with the satellite retrievals being 1 • whilst the individual tower observations are for scales of the order of 1 km or less.The tower R N are from the broadband net radiometer sensors located on each tower.In the absence of reliable measures of G at the tower scale, and in order to derive genuinely independent measures of against which to evaluate the satellite data, we have taken the tower net available energy as the sum of the measured sensible and latent heat flux, i.e.Eq. ( 1).Thereby we have assumed that the eddy covariance flux measurements are able to close the energy balance (i.e.R N − G = λE + H ), the implications of which will be discussed below.We have chosen 30 sites covering a broad range of geographical locations selected from 7 land cover types, including evergreen broadleaf forest (EBF), mixed forest (MF), evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), savanna (SAV), grassland (GRA) and cropland (CRO).A comprehensive list of the site characteristics are provided in Table 1.Each tower evaluation data set is comprised of the 13:30 LT time samples of R N , H and λE which correspond to the satellite overpass.Again, the evaluation is based on pooling these data into weighted monthly average values.For the evaluation we have elected to compare all 12 months of data for 2003, as this year had the best overlap between the FLUXNET and AIRS databases.However before directly validating the satellite-retrieved , the proposed retrieval method is first evaluated using high-temporal-frequency ground-based observations of R N and T S over some eddy covariance sites representing four broad biome categories (grassland, cropland, forest and savanna).Both R N and T S at 13:30 LT and 01:30 LT were extracted from half-hourly measurements, and at 13:30 LT was determined using Eq.(10a, b).The retrieved midday was validated against tower-observed latent and sensible heat fluxes.

Results
Table 2 shows the results from the sensitivity analysis.For G we see the importance of the longwave specification and in particular ε C .The standard deviation of the estimate of G from the ensemble is 18 W m −2 , giving approximate 95 % confidence detection limits of ±36 W m −2 on the estimates.R N is sensitive to all components of the radiation balance calculation as expected (Table 2).The standard deviation of the estimate of R N from the ensemble is 40 W m −2 , giving approximate 95 % detection limits of ±80 W m −2 on the estimates (Table 2).Not surprisingly, the sensitivity results for  mirror those of R N , albeit with a marginally higher ensemble standard deviation of 44 W m −2 (Table 2).
The locations of the 30 terrestrial evaluation sites are marked in Fig. 1. Figure 5 shows annual average, global satellite scenes for 13:30 LT R N , c, G and for the year 2003.Missing data in the images are mainly due to missing data in the AIRS soundings at high latitudes or over the mountain belts, where it is difficult to profile air temperature and relative humidity reliably.In addition, persistent cloudy conditions also prevent reliable retrieval and hence are rejected, although these will be less evident in the monthly or annual average data.
Figure 5a shows the global distribution of R N , which generally decreases with latitude, as expected.R N also decreases over land due to the generally higher albedo, resulting in reduced absorption of the net shortwave radiation (Giambelluca et al., 1997;Gao and Wu, 2014) or relatively higher surface temperature increasing the net longwave component, especially over the drier regions (Liang et al., 1998;Trenberth, 2011).As a result the magnitude of R N was around 200-300 W m −2 over the dry desert regions, whereas the oceanic values of R N were 450-700 W m −2 .
Figure 5b shows the global distribution of c (surface heat capacity).The oceanic values of 4 to 8 MJ m −2 K −1 are equivalent to 1 to 2 m of seawater, which appears reasonable on the daily time step to which they relate (Stramma et al., 1986;Schwartz, 2007).These oceanic values are somewhat noisy due to the small day-night temperature differences observed for the oceans giving a relatively poor signalto-noise ratio.However, behind this noise the pattern of oceanic c appears relatively uniform as one might expect.higher values than the drier, less vegetated areas, as expected.
The soil equivalent depth of this heat capacity is approximately 0.01 m, which again appears reasonable for a daily time step (Li and Islam, 1999), although in heavily vegetated areas c is obviously comprised of a more complex aggregation.
Figure 5c shows the global distribution of G.These are the 13:30 LT values, hence their being net positive as an annual average.Between 20 • north and south, G is approximately 10 to 20 % of R N , and this rises to more than 40 % above 50 • north and south (Hsieh et al., 2009).Given that this opposes the pattern of R N , one would conclude that there are either some deficiencies in the way G is specified here or that R N partitions into latent heat far more effectively than surface heating in these warm wet environments (Liu et al., 2005).Again, terrestrial values are lower than their oceanic equivalents, mainly due to the lower heat capacity as well as reduced R N as discussed above.This also highlights the role of the vegetation layer in preventing ground heating (Baker and Baker, 2002;Bounoua et al., 2010).The Sahara appears particularly prominent in this scene, with high rates of midday heat accumulation, which appears to be associated with a combination of moderate net radiation and relatively high heat capacity.The heterogeneity in this region appears to be related to the pattern of bare darker rock.
Figure 5d shows the global distribution of , which follows a similar pattern to R N as expected, although the pattern of G shown in Fig. 5c dictates that the north-south gradients in are somewhat stronger than those of R N .Before discussing these results, we consider their evaluation.In the first step, we validated the new method of retrieval at representative FLUXNET sites using ground observations of the surface radiation components (T S , R N and G) as input before directly evaluating the satellite-based retrievals.Tower-scale evaluation of daily midday (13:30 LT) is illustrated in Fig. 6a, b, c and d    to 0.98 (±0.04)] 1 between observed and predicted across all the biomes with regression statistics ranging between 0.89 (±0.07) and 1.12 (±0.05) for the gain and −44.29 (±20.15) to 59.40 (±29.03) for the offset (Fig. 6).The root-meansquare deviation (RMSD) varied between 41 (savanna) and 88 W m −2 (forest).
Figure 7a shows the pooled evaluation of satellite R N , which produced an overall correlation of r = 0.88 (±0.03).Assuming both tower and satellite observations are linearly related through some "true" value, then the pooled values are co-related by R N (satellite) = 0.75(±0.02)RN (tower) + 23.37(±8.20),i.e. a small but significant underestimation in R N (satellite) relative to R N (tower).The RMSD between the two was 98 W m −2 .The biome-specific statistics for R N are given in Table 3, which reveals correlations ranging between 0.65 (EBF) and 0.96 (ENF), RMSD ranging between 74 (GRA) and 127 W m −2 (EBF), and regression statistics ranging between 0.58 (±0.08) and 0.87 (±0.04) for the gain and −32.40 (±23.73) and 107.45 (±39.93) for the offset.
Figure 7b shows the evaluation for , which produced pooled statistics of r = 0.87 (±0.03) and an RMSD of 72 W m −2 , and the regression between the satellitepredicted and tower-observed produced a regression line of (satellite) = 0.90(±0.03)(tower)−2.43(±8.19).The biome-specific statistics for are also given in Table 3 showing correlations ranging from 0.70 (EBF) to 0.95 (ENF), RMSD ranging between 62 (GRA & SAV) and 88 W m −2 (EBF) and regression coefficients ranging between 0.66 1 All uncertainties are expressed as ± 1 standard deviation unless otherwise stated.For details of the site characteristics see Table 1.For the comparative statistics see Table 3.The solid line is the pooled linear regression given in Figure 8 shows a sample of monthly time series for for both the satellite and the towers.The sites were selected to represent the biome classes considered here, as well as ones for which complete annual data sets for 2003 were available.These results show that the satellite estimates generally track the trends in the tower data and hence that the pooled statistics are not masking the within site variability.Again, the site-wise comparative statistics for these data are given in Table 3.

Discussion
For R N the statistics relating the satellite and tower data are different with to results of the following authors: Bisht et al. (2005), who obtained 74 W m −2 RMSD when evaluating MODIS Terra geophysical land products over the Southern Great Plains of the USA (our RMSD in grassland is only comparable here, while other biomes show larger error); Jacobs et al. (2004), who obtained a 14-46 W m −2 RMSD (12.2 % relative RMSD) when determining hourly R N using GOES (Geostationary Operational Environmental Satellite) data over wetlands in southern Florida; Cai et al. (2007), who obtained 13.7 % error when evaluating MODIS Terra-Aqua data over China; Bisht and Bras (2010), who obtained 39-51 W m −2 RMSD over the central USA using MODIS Terra atmospheric data at 5-10 km spatial resolution; and Hwang et al. (2013) and Hou et al. (2014), who reported RMSD in instantaneous and daily R N to be 58-142 and 37-40 W m −2 over Southeast Asia and China, respectively, using MODIS Terra data products.Stackhouse et al. (2000)  data (a comprehensive list of relevant studies is given in Table 4).It is important to emphasise that, in the present study, the RMSD is being impacted in two ways: due to spatial scale mismatch and due to the time integration.When these errors are compounded in the derivation of R N and compared with tower data, an RMSD of the order of 98 W m −2 appears reasonable considering the coarse spatial resolution of the AIRS data (1 • × 1 • ).There have been very few attempts to retrieve satellite estimates of and compare these with ground truth data, although the statistics from our attempt appear to be parallel to the results of Stisen et al. (2008), who studied a single grassland site in the Senegal River basin using moderate (high) spatio-temporal resolution (5 km spatial resolution, 15 min temporal resolution) MSG (Meteosat Second Generation) geostationary satellite data and obtained a correlation of r = 0.71 and an RMSD of 43 W m −2 in comparison to the surface measurements.While estimating evapotranspiration over Indian agroecosystems, Bhattacharya et al. (2010) obtained an RMSD of 56 W m −2 for noontime using 8 km resolution Indian geostationary satellite data.In another study with MODIS Aqua data over semi-arid agroecosystems in India, Bhattacharya et al. (2011) reported an RMSD of 34 W m −2 in daily average , which was associated with a significant tendency to underestimate .As was common to all these studies, the ground heat flux was either modelled as an empirical approximation employing remotely sensed surface variables (albedo, vegetation index and T S ) or as a fixed fraction of R N .Murray and Verhoef (2007) argued that these empirical approaches do not generalise well.In particular, prescribing G as a fixed fraction of R N overlooks the role played by the thermal inertia of the land surface (Santanello and Friedl, 2003), leading to an underestimation of G in the morning and overestimation during the afternoon (Gentine et al., 2007).The retrieval of G proposed here using day-night surface temperature information attempts to account for this thermal inertia effect, and the results appear to support this approach especially when considering the scale mismatch between the tower and satellite observations.
As seen in Fig. 7a, there is a systematic underestimation of R N relative to the tower values which exceeds the typical accuracy of net radiometer measurements of 20 W m −2 quoted by Foken (2008).We examined this underestimation in more detail wherever possible by evaluating three of the individual radiation components of R N (R S↓ , R L↓ and R L↑ ).All tower sites provided measurements of R S↓ (but not R S↑ ). Figure 9a shows R S↓ is systematically underestimated at the satellite scale, with R S↓ (satellite) = 0.70(±0.02)RS↓ (tower) + 68(±12.24),which accounts for the mismatch of R N (satellite) ≈ 0.75R N (tower).Before attempting to account for the various reasons for this underestimation, it is important to realise that, unlike the IR components, the shortwave components are all-sky retrievals, i.e. like the tower data they do not omit cloudy-sky conditions.As a result, any bias in the shortwave is not as a result of biased sampling when comprising the monthly average.Furthermore, the omission of non-clear-sky data would tend to lead to R S↓ (satellite) > R S↓ (tower).
Clearly, the retrieval of atmospheric shortwave transmissivity (τ A ) using cloud cover fraction is the principal reason for R S↓ (satellite) < R S↓ (tower) (Fig. 9a).The sensitivity analysis presented in Table 2 also indicates the significant sensitivity of R N and to the cloud cover fraction and atmospheric transmissivity.This shows that the method presented in the manuscript to estimate R S↓ needs further improvements.If we assume τ A to be the principal reason for R S↓ (satellite) < R S↓ (tower), then a global value of 0.75 would be, on average, too low (Gueymard, 2003).A recent study of Longman et al. (2012) for the Mauna Loa Observatory (MLO) demonstrated the clear-sky τ A could go up to 0.90.Given the relatively well defined relationship between R S↓ (satellite) and R S↓ (tower) seen in Fig. 3c, one would imagine that a more sophisticated dynamic representation of τ A would offer substantial improvements in R S↓ (satellite).Retrieval of τ A including other atmospheric (e.g.cloud optical depth, aerosol optical depth, total precipitable water etc.) and surface (for example, single scattering albedo) variables in addition to the cloud cover fraction would offer a potential possibility of refining the R S↓ estimates (Chen et al., 2014;Longman et al., 2012    Hogue, 2008).The exo-atmospheric shortwave radiation frequently interacts with the clouds, aerosols and water vapour during the transmission towards the Earth's surface.This interaction is wavelength-dependent over the entire shortwave spectrum (Chen et al., 2014;Kim and Hogue, 2008;Gueymard, 2003), and therefore a spectrally resolved τ A scheme will be valuable to accurately determine R S↓ .Recalibration of τ A using the tower data is also a possibility, although we have avoided this given that the AIRS cloud cover fraction and scale mismatch between the satellite and tower could also be involved in the observed bias.For example, the diffuse fraction of R S↓ (tower) can become enriched by surfacereflected solar radiation, particularly in undulating terrain (Dubayah and Loechel, 1997;Sultan et al., 2014).Nonlin-ear scaling effects of surface albedo (Oliphant et al., 2003;Salomon et al., 2006) can also be included in this because surface albedo interacts nonlinearly with surface characteristics such as surface wetness and land surface temperature (Ryu et al., 2008) or the leaf area index (Hammerle et al., 2008).Although, the RMSD of instantaneous R S↓ obtained in the present study (110 W m −2 ) is different to other studies where R S↓ retrieval was based on either using parametric (radiative transfer) models or look-up tables derived from high-spatial-resolution MODIS data, it is worth comparing it with the statistics of some of those studies.To probe the specification of R N further, we investigated the individual longwave radiation components in relation to measures of theses fluxes available for a limited subset ( 14) of tower sites where the longwave radiative flux components were directly measured by pyrgeometers.From Fig. 9b and c it appears that there is quite good agreement between the satellite and tower data for both R L↓ and R L↑ and that any mismatch is insufficient to explain the discrepancy in R N .This is somewhat surprising for two reasons.Firstly, unlike the shortwave component, R L (tower) is all-sky whilst R L (satellite) is only from clear-sky conditions where IR retrieval is possible.As a result, one would anticipate very significant differences in the monthly average values of the longwave components.However, it is difficult to predict the effect of this biased sampling on R NL (satellite) given that cloud interacts with both R L↓ and R L↑ in complex ways.Secondly, one would anticipate significant scaling effects from the T 4 nonlinearity in Eq. ( 6), which can result in a disproportionate contribution of warmer elements within the system to both R L↓ and R L↑ (Kustas and Norman, 2000;Lakshmi andZehrfuhs, 2002, Corbari et al., 2010).The fact that these effects are not seen to any significant degree could point to compensating errors in the analysis but does not distract from the central message of the importance of the bias in the shortwave when accounting for R N (satellite) < R N (tower).
Figure 7b and Table 3 show that (satellite) ≈ 0.90 (tower), suggesting a slight compensation for the underspecification of R S↓ through the underspecification of G from the satellite data.However, this evaluation assumes the energy balance to be closed in the tower data (i.e.R N − G = λE + H ), which typically is not the case, with λE + H often falling short of R N − G by 20 % (Wilson et al., 2002).Because the causes of this energy imbalance remain controversial (see Foken, 2008, for a review), it is difficult to estimate how much the tower values of λE + H are actually biased and hence the extent to which this bias affects our evaluation.Stoy et al. (2013) recently found a systematic relationship of the surface energy balance closure with landscape heterogeneity over 173 FLUXNET tower sites and reported an energy imbalance of 9-30 % over diverse biomes.In another study, Amiro et al. (2009) found relatively better fulfilment of energy balance closure by averaging data over longer periods.The monthly averages of (AIRS overpass time) 13:30 LT surface energy balance closure of the 30 sites used here (Table 5) show an average energy imbalance of ∼ 20 % (ranging from 8 to 34 %).The errors in the tower data are also believed to be associated with different footprint characteristics for the different instruments used (Lin et al., 2008).For example, R N observations typically have a footprint size of ∼ 10 m 2 , whilst air properties (e.g.air temperature, humidity) have footprint sizes of > 1 km 2 .By way of illustration, if, in the worst case, the entire energy imbalance were to be attributed exclusively to λE + H (i.e.R N − G are quantified correctly), then the true midday λE + H could be some 20 % greater (Wilson et al., 2002).As a result, the present bias seen in Table 3 would change to (satellite) ≈ 0.72 (tower) , again suggesting R S↓ to be the main source of bias in the satellite retrievals for both R N and .
We have only evaluated the satellite retrievals using data from terrestrial sites, and clearly it would be worthwhile repeating this for the ocean retrievals if possible.We have held back on this evaluation here because of the lack of an extensive network of instantaneous latent and sensible heat flux or radiative flux data over the oceans, although we note that the SEAFLUX project within the Global Energy and Water Experiment (GEWEX) initiative should give rise to such a database in the near future.From the terrestrial evaluation we would argue that the methodology employed here shows promise for specifying both noontime R N and , although the results suggest the need for improvements, particularly in the specification of R S↓ .More detailed studies evaluating the representativeness of each tower site footprint in relation to the 1 • scale within which it is situated could prove useful in this regard as would methods for cloud-proofing the satellite retrievals under persistent cloudy-sky conditions.Similarly, evaluation under extreme conditions (e.g.high altitudes and latitudes) is also required.

Conclusions
We have demonstrated a novel retrieval for midday (13:30 LT) surface net available energy ( ) by blending the monthly atmospheric and land surface variables of AIRS and MODIS sensors.We have attempted to structure the method such that the retrieval does not depend on any offline calibration.We performed a two-step evaluation of the retrieved values for at some representative surface radiation measurement sites and also over 30 FLUXNET sites from all over the globe.
Current estimates performed well when compared against high spatio-temporal in situ measurements and coarse-spatial-resolution satellite retrievals.Consistent underestimation of R N was noted due to the underestimation of shortwave radiation, and in tropical latitudes the R N (and ) agreement was relatively weaker.Combination of highfrequency cloud dynamics and relatively low seasonal variability of R S↓ makes it difficult to accurately model both and R N .One of the key challenges in modelling R N (and ) in the tropics is to account for fast-changing cloud cover fraction and cloud optical properties throughout the day.
With the availability of high-spatial-resolution (1-5 km) MODIS day-night optical and thermal data, our present approach could be extended to derive high-spatial-resolution estimates at the global scale.This could be achieved by estimating MODIS-based day-night R N and combining day-night R N with day-night T S observations.At the same time, the current methodology could also be used on hightemporal-frequency observations of geostationary satellites (e.g.GOES and METEOSAT).Once surface heat uptake and heat capacity (through Eq. 10) has been estimated, hourly G and could be determined from hourly R N and T S observations of geostationary satellites by assuming conservation of heat capacity over a particular day.Operational generation of satellite-based product would be a valuable resource for a variety of investigations, such as estimating latent and sensible heat, evaluation of Earth system model outputs and quantifying the land-atmosphere coupling strength.Given that we have resorted to the minimal amount of calibration in deriving , it would appear sensible to adopt a similar philosophy in developing satellite-based schemes for latent and sensible heat fluxes as proposed in M2.
In addition to opportunities to specify large-scale surface heat and water vapour fluxes, the heat capacity estimates made here clearly carry information on variations in terrestrial properties such as surface moisture storage, and we envisage that studies to develop this concept further could prove fruitful, particularly because of the emergence of satellite microwave data against which the results could be compared.

Figure 1 .
Figure 1.The distribution of the 30 eddy covariance tower sites used for evaluating R N and .

Figure 2 .
Figure 2. Illustrative examples of in situ monthly diurnal surface temperature (T S ) and the saw-tooth pattern of monthly satellite midday (13:30 LT) to night-time (01:30 LT) T S (as hypothesised in Eq. 10) for three different seasons of a year.This shows a linear rise and fall of day-night T S (dotted black line) or vice versa and indicates T S symmetry.This also shows how well the two T S samples capture the dynamic range of the day and hence the discretisation is representative of the daily energy balance.The examples in Fig. 2 are shown with half-hourly data.

Figure 3 .
Figure 3. Examples of the diurnal evolution of net available energy ( ) (= λE + H ) during three different times of year over four representative biome types.This shows ∼ = 0 around 01:30 LT.
Over land c varies between 0.05 and 0.5 MJ m −2 K −1 with wetter tropical and high-latitude areas showing significantly

Figure 7 .
Figure 7.Comparison of satellite and tower monthly average 13:30 LT (a) R N and (b) [ = (H + λE) (for the tower sites)].For details of the site characteristics see Table1.For the comparative statistics see Table3.The solid line is the pooled linear regression given in Table3.(+ EBF; × MF; • GRA; * CRO; ∇ ENF; DBF; SAV.) evaluated the International Satellite Cloud Climatology Project (ISCCP) data and found errors in the range of 10 to 15 W m −2 in monthly average shortwave and longwave radiative fluxes.Other studies reported 33-60 W m −2 RMSD in daily R N using 5 km MODIS Terra optical and thermal www

Figure 8 .
Figure 8. Satellite (grey) and tower (black) time series of monthly average 13:30 LT net available energy [ = (H + λE) (for the tower sites)] for a selection of sites for 2003.The numbers on the x axis are the month numbers, i.e.January is month number 1 and December is month number 12.

Table 1 .
Eddy covariance sites used for the evaluation of the satellite-derived R N and .

Table 2 .
Sensitivity analysis results.The forcing data are taken for midsummer in the Southern Great Plains, USA.Sensitivities are locally linear, averaged across the ensemble response and expressed as dimensionless relative changes.Only absolute sensitivities > 0.1 are shown.N = 10 5 realisations.

Table 3 .
Comparative statistics for the satellite-and tower-derived monthly midday (13:30 LT) R N and for a range of biomes.Values in parenthesis are ± 1 standard deviation.

Table 4 .
Summary of errors and characteristics of some of the dedicated satellite-based R N and R S↓ retrieval studies.

Table 5 .
Stoy et al. (2013)ard deviation of monthly midday (13:30 LT) surface energy balance closure (C EB ) of different biome types for the 21 (out of 30) FLUXNET sites (N) under study.Nine out of 30 sites had missing ground heat flux.C EB is computed according toStoy et al. (2013).−2 at 1 • spatial resolution for the instantaneous R S↓ estimates and 20-39 W m −2 for the daily R S↓ (and net shortwave, R NS ) estimates.Considering the simplicity of the current approach and the large spatial scale of the AIRS data, an RMSD of the order of 110 W m −2 appears reasonable.