Combining satellite observations and reanalysis energy transports to estimate global net surface energy fluxes 1985-2012

14 Two methods are developed to estimate net surface energy fluxes based upon satellite-based 15 reconstructions of radiative fluxes at the top of atmosphere and the atmospheric energy tendencies 16 and transports from the ERA-Interim reanalysis. Method 1 applies the mass adjusted energy 17 divergence from ERA-Interim while method 2 estimates energy divergence based upon the net 18 energy difference at the top of atmosphere and the surface from ERA-Interim. To optimise the 19 surface flux and its variability over ocean, the divergences over land are constrained to match the 20 monthly area mean surface net energy flux variability derived from a simple relationship between the 21 surface net energy flux and the surface temperature change. The energy divergences over the oceans 22 are then adjusted to remove an unphysical residual global mean atmospheric energy divergence. The 23 estimated net surface energy fluxes are compared with other data sets from reanalysis and 24 atmospheric model simulations. The spatial correlation coefficients of multi-annual means between 25 the estimations made here and other data sets are all around 0.9. There are good agreements in area 26 mean anomaly variability over the global ocean, but discrepancies in the trend over the eastern 27 Pacific are apparent.


Introduction
The absolute mean value of net radiation imbalance at the top of atmosphere (TOA) is a key climate variable, providing an estimate of total energy gain of the Earth system and a link between radiative forcing, ocean heat uptake and surface temperature response.It has been estimated to be 0.5 to 1 W/m 2 for the global mean in recent studies [Hansen et al., 2011;Loeb et al., 2012;Trenberth et al., 2014;Wild et al., 2015] using changes in total ocean heat content (OHC) [Lyman and Johnson, 2014;Trenberth et al., 2014;Smith et al., 2015;Roemmich et al., 2015] and making assumptions about minor energy sinks including the land, the atmosphere and the cryosphere.Although satellite data provide regional coverage of top of atmosphere radiative fluxes, the net surface fluxes display much larger uncertainty due to the lack of constraints from global observations [Trenberth et al., 2009;Wild et al., 2013;Wild et al., 2015].
The net energy fluxes at the earth's surface, including short and long-wave radiation and the sensible and latent heat fluxes, are very important for the study of surface temperature change, and the atmospheric and oceanic circulations.The surface fluxes also control the water cycle since the incoming short wave radiation provides much of the energy required for surface water evaporation.
Net downward surface energy can accumulate within the ocean, leading to long time-scale effects on the climate.Therefore accurate estimation of the surface energy fluxes is essential for understanding both the short term temperature hiatus [Easterling and Werner, 2009;Knight et al., 2009;Trenberth and Fasullo, 2013a;Huber and Knutti, 2014;Watanabe et al., 2014] and long term climate change [Otto et al., 2013].It is difficult to obtain accurate absolute surface fluxes from satellites due to complicated atmospheric conditions affecting the retrieval processes in particular relating to the numerous surface variables required by turbulent flux bulk formulae [Schmetz, 1991].
The net input of radiation fluxes at TOA are modulated by the atmosphere and re-distributed by lateral energy transports [Keith, 1995;Chiodo and Haimberger, 2010;Lucarini and Ragone, 2011;Trenberth and Fasullo, 2013a;Mayer and Haimberger, 2012;Mayer et al., 2014;England et al., 2014;Loeb et al., 2015].Meehl et al. [2011] and Trenberth and Fasullo [2013b] also demonstrated that the vertical energy redistribution in the oceans is likely to have contributed substantially to the slowing in the rate of global average surface temperature increase in the last fifteen years.
Assessment of where the net accumulation of energy in the climate system is being stored within ocean basins [Balmaseda et al., 2013;Drijfhout et al., 2014;Llovel et al., 2014;Desbruyères et al., 2014;Roemmich et al., 2015] is required for understanding the mechanisms of energy redistribution associated with internal variability and therefore the surface temperature variations.
The currently available surface flux data sets have some limitations.Observed data from in situ measurements are sparsely distributed in space, while satellite-derived retrievals contain substantial uncertainties and require further validation.Observationally-based data, reanalysis estimates and climate model simulations show a large spread in the data and large unrealistic global imbalances when turbulent and radiative flux products are combined [Trenberth et al., 2009;Stephens et al., 2012;Wild et al., 2013].In this study, we apply an atmospheric energy divergence approach [Chiodo and Haimberger, 2010;Mayer and Haimberger, 2012] using two different methods to estimate the net downward surface energy fluxes by combining reconstructed net radiation fluxes at TOA [Allan et al., 2014] with the energy tendencies and lateral divergence simulated by the ERA-Interim reanalysis [Dee et al., 2011;Berrisford et al., 2011].

Data sets
The key data set is the ECMWF (European Centre for Medium-Range Weather Forecasts) ERA-Interim reanalysis (ERAINT) [Dee et al., 2011].Various observational data are assimilated to a weather forecast model to provide representations of atmospheric states.Although it has some known problems, such as the lack of volcanic aerosols and the omission of the 11 year solar cycle [Dee et al., 2011], it provides a comprehensive representation of atmospheric variables and estimates of energy divergences and fluxes required for this study.The net radiation flux at TOA is based on the recent reconstruction by Allan et al. [2014] using satellite observations from the Clouds and the Earth's Radiant Energy System (CERES; Loeb et al., 2012) and Earth Radiation Budget Satellite (ERBS) wide field of view (WFOV, 72 day mean) non-scanning instrument [Wong et al., 2006], ERA-Interim reanalysis and climate model simulations applying the Atmospheric Modelling Intercomparison Project 5 (AMIP5) experimental setup with prescribed observed sea surface temperature (SST) and sea ice and realistic radiation forcings [Taylor et al., 2012].The net TOA flux is adjusted to ensure agreement with an observational estimate over the period 2005-2010, primarily determined by observed 0-2000m ocean heating rate [Loeb et al., 2012;Allan et al., 2014].The TOA reconstructions are updated using the latest version (version 2.8) of CERES data.Another important update from Allan et al. [2014] is that prior to March 2000, reconstructed radiative fluxes are adjusted separately for each hemisphere rather than applying a global adjustment.This adjustment ensures that deseasonalized anomalies in radiative fluxes match the WFOV variability for 0-60 o S and 0-60 o N regions.Further details of the additional adjustment procedures are described in Allan et al. [2014].The updated net downward TOA radiation flux will be referenced as   .
Sixteen AMIP5 models are used in this study and one member from each model is chosen.Data from a 5 member ensemble of the UPSCALE (UK on PRACE -weather-resolving Simulations of Climate for globAL Environmental risk) [Mizielinski et al., 2014] simulations are also used here.UPSCALE is from a global atmospheric model (HadGEM3-A-GA3; Walters et al. [2011]) at 25km resolution, which is employed to produce an extended AMIP simulation up to 2011 using the Operational Sea Surface Temperature and Sea Ice daily high resolution Analysis (OSTIA, Donlon et al. [2012]).The only differences between these 5 ensemble member runs are their initial conditions: each member was perturbed by randomly altering the lowest order bit in the 3D potential temperature field.
The recently available ECMWF 20 th century atmospheric reanalysis from ERA-CLIM (European Reanalysis of Global Climate Observations) project (hereafter ERA20C) is used here for comparison purpose; it is a single member reanalysis and it assimilates observations of surface pressure and surface marine winds; SST, sea ice and realistic radiative forcings are prescribed [Poli et al., 2013].
The atmospheric energy divergence from the MERRA (Modern Era-Retrospective Analysis for Research and Applications) reanalysis is also used for the comparison of net surface energy fluxes.A large quantity of observational data are assimilated in the MERRA system using a three-dimensional variational data assimilation analysis algorithm [Rienecker et al., 2011].Observed surface temperature data are from HadCRUT4 [Morice et al., 2012].All data used in this study are monthly mean diagnostics accumulated from higher time resolution data and are listed in Table 1.

Surface energy flux from mass adjusted divergence
Following Berrisford et al. [2011], the total energy (E) in an atmospheric column can be written as where L, q, Cp, T,   and k are the latent heat of condensation of water, specific humidity, the specific heat capacity of air at constant pressure, temperature, surface geopotential and kinetic energy (( • )/;  is the horizontal wind velocity vector), respectively. is the pressure,  is the gravitational acceleration and η is the hybrid vertical coordinate which is a function of pressure and surface pressure [Simmons and Burridge, 1981].The total energy tendency is the energy divergence (  ).The horizontal flux in   is not simply the flux of total energy from equation ( 1), but incorporates the flux of enthalpy [Boer, 1982;Trenberth and Solomon, 1994).
For mass consistency, the output   from ERA-Interim should be mass adjusted, because during the assimilation procedure, observations reset the surface pressure field, whereas the mass fluxes are not adjusted accordingly [Graversen et al., 2007;Berrisford et al. 2011 , where   is vertically integrated water vapour divergence and   is total column water vapour tendency which can be calculated from the time series of total column water vapour content.Both are obtained from ERA Interim; this method is considered more accurate than using   −  directly from the reanalysis, since water vapor is assimilated, but precipitation is a simulated variable that is highly dependent upon model parameterisations.It includes water mass transfer due to phase change between water vapour and liquid water.The phase change between liquid water and ice in the atmosphere has been ignored and the horizontal water transport due to cloud advection is also neglected since these terms are small, ℎ ̅ and  ̅ are the vertical average of moist static energy and kinetic energy, respectively, which can be computed from analysed ERA-Interim fields.
From equation (3), we can have the net surface energy flux from the mass adjusted energy divergence (  ): Similar procedures are applied to MERRA data [Mayer et al., 2013] to obtain mass adjusted total energy divergence  − which is substituted into equation ( 6) to obtain the net downward surface flux  − .

Surface energy flux from model residual divergence
Another way to estimate the atmosphere energy divergence is to calculate it directly from ERA-Interim as a residual of energy fluxes [Chiodo and Haimberger, 2010;Mayer and Haimberger, 2012]: , the effect of aerosolrelated biases on   will be reduced.

Adjustment constraints
Since a large quantity of observational data are assimilated into ERA-Interim, it is expected both   and   will provide reasonable spatial structures, but the   has a multi-annual (2001)(2002)(2003)(2004)(2005) global mean value of -0.9 Wm -2 which is not physically reasonable since it is expected the global averaged   should be zero to guarantee energy conservation.This is because atmospheric models don't, in general, have a closed budget for the atmospheric energy, as a result of inconsistent treatment of turbulent cascades of kinetic energy and water mass [Lucarini and Ragone, 2011;Previdi and Liepert, 2012;Lucarini et al., 2014].over land and ocean respectively and there is a net energy transport from the ocean column to land column [Wild et al., 2015].The steps for estimating the monthly net surface energy fluxes are as follows: Remove the global mean divergence from   .
We already have   and   , assuming we have the correct monthly net surface energy flux data over land, the monthly vertically integrated energy divergence can be calculated over land using energy balance equation; The globe is divided into 15 o latitude band (30 o over Antarctic).The mean discrepancy between mass corrected divergence and the one derived from step (2) over land is redistributed evenly over ocean grid points to keep the total divergence unchanged across each band.
The monthly net surface energy flux over the ocean can then be calculated using bias corrected divergence.
In step (2), it will be ideal to use net surface energy flux calculated from   as the initial estimation over land, but as mentioned above the derived fluxes have unrealistically large regional changes (2001-2008 mean minus 1986-2000 mean) over land, so the surface energy flux from ERA-Interim ( − ) over land is used as the initial estimation.In order to correct the unrealistic trend and large anomaly variability of  − as discussed in section 3.3 (which would imply large unrealistic temperature variations or land heat capacity), a simple method described in the next section is applied to estimate the monthly net energy flux variability based on the observationally constrained surface temperature changes over land.

Net energy flux over land
The mean global land flux is estimated using the simple relationship of   =  Based upon the 1985-2012 mean surface temperature change of 0.0298K/year from HadCRUT4 we estimate the mean of the reconstructed net surface flux as 0.08W/m 2 over this period.Setting this flux to zero is also reasonable [Trenberth et al., 2009].Combining algorithms in sections 2. The difference is mainly due to the albedo difference between the land and the ocean.The large energy deficit over land should be compensated by the horizontal energy transport from ocean to land [Mayer and Haimberger, 2012;Trenberth and Fasullo, 2013b].

Net energy flux at the surface
The multi-annual mean (2001)(2002)(2003)(2004)(2005) net surface energy fluxes from   are plotted in Fig. 3a and zonal mean variations from   ,   ,  − , ERAINT, ERA20C, UPSCALE and AMIP5 data sets are plotted in Fig. 3b-d.The area-weighted means are displayed in the zonal mean plot.The multi-annual mean for other data sets are in Fig. S2.  and  − are calculated from the spatially filtered   and  − respectively using a Hoskins spectral filter [Sardeshmukh and Hoskins, 1984] with an attenuation of 0.1 at wave number 106 [Berrisford et al., 2011].A filter is necessary due to the noise generated by data assimilation, highlighting that spatial patterns must be interpreted with caution.Ragone [2011] and coincide with the maximum in baroclinic activity [Lucarini and Ragone, 2011].
The transport from  − has stronger magnitude at 40 o S/N compared with the other estimates.The transport from   is of larger magnitude than that from   in the northern and southern hemisphere sub tropics, consistent with Mayer and Haimberger [2012].
Due to flux constraints over land, the area mean fluxes from both   and   are identical.Their spatial structures and zonal mean variations are also very similar (Fig. 3 and Fig. S2a), but the magnitudes differ in places as shown in Fig. 5a.  is larger in magnitude than   in the south Indian Ocean, but smaller in the north Indian Ocean.  is smaller over the central, west and north west Pacific, but has larger values over the subtropical gyre of north Pacific, as well as over south east Pacific.
Though the mean surface flux spatial structure of ERAINT (Fig. S2b) is similar to the derived ones, its area mean fluxes are unrealistically large over the global ocean (9.30Wm -2 in Fig. 3c) compared with ocean observations [Llovel et al., 2014;Roemmich et al., 2015] which are of the order of 0-1 W/m 2 .ERA-Interim surface fluxes are substantially larger than   over the oceans as shown in Fig. 5b, except for the area near the Equator, and this can be seen clearly from the zonal mean variations (Fig. 3c).ERA20C simulates large fluxes into the Southern Ocean, more flux from ocean to atmosphere over the whole Indian Ocean and the north and south Atlantic subtropical gyres (Fig. 5c).As stated earlier, the  − (Fig. 5d are not apparent in other data sets (Fig. 5e).The ensemble mean from AMIP5 simulations show much lower fluxes into the Western Pacific (Fig. 5f) and this is mainly contributed from CMCC, CNRM, FGOALS, GISS, MRI and INMCM4 model simulations as shown in Fig. S3.

Changes in downward energy flux
In order to investigate where the energy is moving through the climate system [Lucarini and Ragone, 2011;Mayer and Haimberger, 2012;Guemas et al., 2013;Allan et al., 2014;Mayer et al., 2014;Drijfhout et al., 2014] MERRA ( − in Fig. 6f) is even noisier than those from   and   , but it also displays decreased net downward energy flux over the eastern Pacific.This has been identified as an important region in determining aspects of the recent slowing rate of global surface temperature rise [Kosaka and Xie, 2013;Trenberth and Fasullo, 2013a;Meehl et al., 2014].On one hand, the cooling Eastern Pacific will suppress turbulent energy transport from ocean to the atmosphere, so the net downward flux would be increased over this region; on the other hand as demonstrated by England et al. [2014], the cooling is due to the observed pronounced strengthening in Pacific trade winds which are not represented fully by AMIP simulations.The increased winds will cause more evaporation, so more latent heat transports to the atmosphere.Brown et al. [2014]  For the reconstructed surface fluxes (  and   ), the global changes from the 1990s to the 2000s (see table S1) are consistent with Allan et al. [2014], who considered the TOA net imbalance; there is an increase in net downward flux at the surface due to the recovery from Pinatubo [Smith et al., 2015].Consistency with global-mean TOA fluxes is expected since the surface flux estimates are based upon these TOA reconstructions and atmospheric heat capacity is small so cannot uptake a significant fraction of the top of atmosphere imbalance [Palmer and MacNeal, 2014].The ocean heat uptake is also increasing since over 90% of the excess energy into the Earth system is stored in the ocean [Trenberth and Fasullo, 2013a].Consistency between global mean surface and TOA flux changes also applies to ERA20C reanalysis, UPSCALE and AMIP5 simulations (see table S1).Smith discontinuities on these changes [Mayer et al., 2013] in the reanalysis, due to artifacts of the observing system, merits further investigation, though the effect is most significant for the partition of the latent and dry static energy and less prominent when considering the total transport [Trenberth and Fasullo, 2013b].
The deseasonalised anomaly (calculated relative to the 2001-2005 period) time series of the area weighted mean net downward energy fluxes at the surface from different data sets are plotted in Fig. 7 for the globe, the global ocean and the global land.The time series from both derivations (  ,   and  − ) are identical by design.The light grey shadings are ± 1 standard deviations of the sixteen AMIP5 simulations.All lines are 6 month running means.The ERAINT data are also plotted for reference purpose; spurious trends are explained by latent heat flux changes over the ocean [Chiodo and Haimberger, 2010] and from longwave radiation over the land.There is good agreement between derived fluxes and those from AMIP5, ERA20C and UPSCALE data sets over the globe.The correlation coefficients between derived and AMIP5, ERA20C and UPSCALE are 0.38, 0.52 and 0.47, which are significant based on the two-tailed test using Pearson critical values at the level of 5%.Over the global ocean, the coefficients are 0.33, 0.52 and 0.45, which are also statistically significant.Over land the correlation coefficient between derived and ERA20C is 0.60.
The correlation coefficients between other data sets can be found in Table S2 and the correlation coefficient maps are in Fig. S4.Future work will consider in more detail the variability across individual ocean basins and comparisons with independent datasets [Drijfout et al., 2014;Mayer et al., 2014;Desbruyères et al., 2014;Roemmich et al., 2015] contributing toward understanding of variation in energy flows into the ocean.

Summary
Surface fluxes are a crucial component of the climate system yet global-scale observational estimates are highly uncertain [Wild et al., 2015].To complement the existing set of surface flux datasets, an alternative method is developed.The net downward energy fluxes at the Earth's surface are estimated through the combination of the reconstructed TOA radiation fluxes [Allan et al., 2014] and the atmospheric energy divergences (Fig. 1) which are calculated using two distinct methods: (1) mass adjusted energy divergence computed from ERA-Interim reanalysis [Trenberth, 2001;Mayer and Haimberger, 2012;Berrisford et al., 2011]; (2) the residual from the difference between the energy fluxes at the TOA and the surface from ERA-Interim.
To correct for unrealistic variability in energy fluxes over the land a correction was applied using a simple mean relation between surface flux and surface temperature change in UPSCALE climate model simulations which are strongly dependent upon the model's land surface component, JULES.
By setting the global energy divergence to zero, applying the corrected surface fluxes over land and adjusting atmospheric energy divergence from the ocean to the land accordingly the net surface energy flux over ocean could be derived.Although this method relies upon the gross relationship between surface temperature change rate and energy fluxes from a simulation and other assumptions it was found that the sensitivity of the ocean surface flux changes to the methods applied over land are relatively small compared to the differences amongst datasets.
The accuracy of the resultant surface fluxes relies heavily on the quality of the reanalysis.The current version of ERA-Interim has some known problems including drifts in energy fluxes and deficient radiative forcing changes relating to anthropogenic and natural aerosol, and problems with mass divergence and conservation [Berrisford et al., 2011].All these will affect the quality of our product.The assimilation of various observed fields into the model draws towards an observed atmospheric state, so the aerosol effect on the mass adjusted energy divergence (  ) should be less than the effect on   , but the accuracy of the divergence relies on other factors too, such as model spin-up and large time sampling errors, as discussed by Chiodo and Haimberger [2010].
Different datasets capture the general global patterns of the multi-annual mean net downward surface fluxes despite the contrasting methods involved.The spatial correlation coefficients of multi-annual means (2001)(2002)(2003)(2004)(2005) between the reconstruction and other data sets are all around 0.9.The area mean surface flux anomaly time series shows reasonable agreement with AMIP5 (r=0.33),ERA20C (r=0.52) and UPSCALE (r=0.45)simulated monthly anomalies over the global ocean.
The change between the 1990s and 2000s over the eastern Pacific differs between datasets: while climate model simulated surface fluxes increase over the period [Katsman et al., 2011], the reconstruction indicates a reduced net downward surface flux.The cooling surface supresses the airsea turbulent energy exchange, but the strengthening of the observed trade winds [England et al., 2014] over this area will reduce the net downward energy flux.Feedbacks involving low-altitude cloud and reflected shortwave radiation may also amplify this response [Brown et al., 2012].Since the estimated surface fluxes are strongly dependent upon the ERA Interim as well as the MERRA reanalysis which both have temporal homogeneity issues [Mayer et al., 2013 ], further verification of these products with other data sets from observations, reanalysis and model simulations is required in order to further understand the strengths and weaknesses of the current methodology.
Assessing the degree to which SST patterns are driving or being driven by surface flux changes in this region merits investigation [Mayer et al., 2014;Drijfout et al., 2014;Desbruyres et al., 2014].
More detailed assessments of recent changes in surface energy fluxes entering distinct ocean basins [Mayer et al., 2014;Desbruyres et al., 2014] will contribute toward improved understanding of energy flows and internal variability in the climate system.

Figure captions
),where  − and  − are energy fluxes at the TOA and surface computed directly from the ERA-Interim 12-hourly forecasts, where their radiation components (shortwave and longwave) are calculated from the radiation transfer model based on the atmospheric states. − also includes turbulent fluxes simulated by the reanalysis.The term( and Haimberger, 2012] and is preferred over analysed tendencies to be consistent with forecasted TOA and surface fluxes.The calculated   can be used to estimate the surface flux (  ) using the reconstructed TOA flux and total energy tendency from ERA-The accuracy of this divergence relies on the accuracy of the atmospheric properties, the radiative transfer through the atmosphere and the turbulent energy calculations at the surface.It is known that ERA-Interim does not represent aerosol forcing due to volcanic eruptions, most notably following the Mt.Pinatubo eruption[Allan et al., 2014], which might affect the divergence (  ) accuracy since the radiation fluxes are affected by aerosols.Although the constraint on divergence is poor, hence the need for mass adjustment, data assimilation constrains parameters towards an observed atmospheric state; with the inclusion of analysis increment, Even though the global mean   is close to zero (~10 -4 ), the net surface flux derived from it has unrealistically large local changes (2001-2008 mean minus 1986-2000 mean, not shown here) and the global mean RMS (root mean square) of the multi-annual mean differences (2001-2008 mean minus 1986-2000 mean) is about 8.5 Wm -2 .The area mean   over land is also large (about 2 Wm -2 over 2001-2005).A strategy was required to address these problems.The schematic diagram shown in Fig. 1 illustrates the energy flow terms used in the estimation of net surface energy fluxes.The left and right columns depict the energy flow where C is the effective mean surface land heat capacity,   is the global land mean surface temperature change rate and  is a constant indicating the energy flux penetrating beneath the surface layer.Data from five UPSCALE ensemble members are used for this estimation.The land surface model in UPSCALE simulations is JULES (Joint UK Land Environment Simulator) which has an explicit representation of the surface energy balance for vegetation, capturing the weaker coupling that exists between the canopy and underlying soil [Best et al., 2011].The effective land heat capacity depends on the soil and canopy properties and the soil water content.After testing we found high correlations between energy flux and the rate of surface temperature change if   is calculated from consecutive months, e.g. the climatology of   in April will correlate well with   calculated from the climatology difference between April and March, so the effective land heat capacity C and the constant  are calculated by regression using the climatology of   and climatological   .The anomaly time series from modeled and reconstructed (from C,   , and ) land surface mean fluxes are plotted in Fig. S1.The correlation coefficients (r) between monthly anomalies (reference period 2001-2005) are all above 0.6.The plotted lines are 6 month running means and the inflated reconstructed lines (red) are multiplied by the ratio of the standard deviation between modeled and reconstructed monthly flux anomalies (values in red in the plot).The variability in   is generally well captured although there are exceptions, notably over the Mt.Pinatubo eruption period since the constant seasonal C is used while in reality it should vary under anomalous situations; as discovered byIles and Hegerl [2014], the models underestimate the precipitation over Pinatubo eruption period which affects the soil moisture content, therefore affecting the relations between temperature change and energy fluxes.Another factor affecting the net surface energy flux variability is the snow and ice melting.While there are considerable limitations, this method was applied to ensure that large biases in the variability in   over land did not diminish the realism of diagnosed   over ocean which is the goal of the present study.Five sets of the regression coefficients from five UPSCALE members using the above method are applied to the global land mean surface temperature (skin temperature) rate of change   from ERA-Interim to get five proxies of mean surface flux; the ensemble mean is used as our estimated global land mean surface net energy flux.Based onBeltrami et al. [2002], the mean net energy flux over the continental lithosphere is 0.0391W/m 2 over 1950-2000, where the mean land surface temperature change from HadCRUT4[Morice et al., 2012] is about 0.0138K/year (from regression).
3 and 2.4, the estimated 2D net surface energy fluxes over land maintains the spatial structure of  − , but the monthly area weighted mean values match those from the simple model (  = downward radiation flux anomalies at TOA are updated fromAllan et al. [2014] using the latest version (version 2.8) of CERES data and adjusting pre-CERES variability to match the interannual anomalies from the WFOV instrument for each hemisphere separately rather than using the 60 o S-60 o N near-global mean.The TOA flux anomaly time series are plotted in Fig.2for the global mean, the global ocean and the global land, respectively.The reference period is from 2001-2005, but WFOV has a reference period of 1985-1999 and is adjusted, for clarity, to match the mean   (reconstruction) anomaly over this period.There is good agreement between variability depicted by   and the other data sets over the global ocean and the globe.The correlation coefficients (r) between   and ERAINT, UPSCALE or AMIP5 monthly anomaly time series are 0.63, 0.60, and 0.58 over the global ocean and 0.64, 0.44, and 0.46 over the land, respectively.All these correlations are significant based on the two-tailed test using Pearson critical values at the level of 5%.The degree of freedom of the time series is calculated by first determining the time interval between effectively independent samples[Yang and Tung, 1998] but additionally assuming that periods separated by 12 or more months are independent.Although ERAINT does not represent changes in aerosol emissions, most notably following the Mt.Pinatubo eruption in 1991, the correlation coefficient between   and ERAINT is still the highest.This reflects the realistic monthly variability of atmospheric circulation patterns through the extensive assimilation of conventional and satellite data by ERA-Interim.The area weighted multi-annual mean net downward energy fluxes from   (Fig.2d) over2001-    2005 are 0.51, 8.35 and -19.0W/m 2 for the globe, the global ocean and the global land, respectively.
Despite the contrasting methods and datasets, the multi-annual means for the period 2001-2005 from all data sets show similar spatial structures and zonal means except for the MERRA data which show much stronger fluxes over the central Indian Ocean and central western Pacific.The spatial correlation coefficients of multi-annual means between estimations and other data sets are all around 0.9.Over the oceans, despite ~10-20 W/m 2 differences present in the zonal means (Fig.3c), all datasets capture the positive downward energy flux over the equatorial central and east Pacific areas due to the interaction between the tropical instability waves[Willett et al., 2006] and the equatorial Pacific cold tongue[Martínez-Garcia et al., 2010] controlled by ocean mixing[Moum et al., 2013].The evaporation is less and there is lower outgoing longwave radiation over this cold region compared with surrounding regions.The negative downward fluxes over the Gulf Stream in the North Atlantic and Kuroshio currents in the North Pacific are due to heat and moisture transport from the warm ocean surface to the cold atmosphere above[Kwon et al., 2010].Over the global land, the UPSCALE simulation has a similar large magnitude residual flux (-0.68W/m 2 ) to the ERAINT flux (0.71W/m 2 ) because it does not have a closed energy budget[Lucarini and Ragone, 2011].This is in part because the high resolution version of the UPSCALE simulations used were not re-calibrated using observations since a key aim of this project was to understand the influence of resolution upon mean climate.The unrealistically large magnitude values at around 55 and 65 o S (Fig.3d) are caused by single grid points at the southern tip of South America and northern tip of the Antarctic peninsula requires further investigation.The mean northward total meridional atmospheric energy transport calculated from   ,   and  − are also plotted in Fig.4a.Peak magnitudes of around 5PW (1PW = 10 15 W) close to 40 o S and 40 o N are broadly consistent withMayer and Haimberger [2012] andLucarini and ) has larger values over the central Indian Ocean and central western Pacific, but smaller values over much of the eastern Pacific.UPSCALE shows the common feature of smaller flux over the north Indian Ocean and larger energy flux over the Southern Ocean, but the strong flux over the western Pacific and smaller energy flux over the Eastern Pacific , considering the changes of multi-annual means in the net downward energy fluxes at both TOA and surface are informative.A preliminary assessment of the multi-annual mean changes (2001-2008 mean minus 1986-2000 mean) from reconstruction (  ,   ,   and  − ), UPSCALE and AMIP5 data sets are presented in Fig. 6.As discussed by Allan et al. [2014], all three data sets show decreased TOA net fluxes over the tropical east Pacific (left column of Fig. 6).The magnitudes of the TOA flux changes over oceans are much smaller than those at the surface.At surface, the estimated changes over land areas are small from estimation (  ,   and  − ), but the flux changes over Russia are slightly larger than in the UPSCALE and AMIP5 simulations.Fig. 6b and Fig. 6d show the increasing downward energy flux over the North Pacific and Southern Ocean (increased ocean heat uptake), but negative flux changes over the central Pacific, north Indian ocean and north Atlantic.Although the individual surface flux components are not reconstructed, considering those simulated by ERAINT, the changes appear to be dominated by latent heat fluxes.Comparing with atmospheric model simulations, although both ensemble means from UPSCALE and AMIP5 simulations show decreased fluxes into the central Indian Ocean and north Atlantic (Fig. 6i, l), the big differences are over the Eastern Pacific where simulated increases in downward flux are opposite to the estimations in Fig. 6b,d,f.The estimated surface flux from also showed that the surface cooling over Eastern Pacific will enhance the reflected short wave radiation, therefore reduce the net downward energy flux.The eastern tropical Pacific region marked in Fig. 6b,d,f covers 20 o N-20 o S and 210 o E to the west coast of the central America.The mean TOA flux change(2001-2008 mean minus 1986-2000 mean)      over this area (Fig.6a) is -2.1W/m 2 while the surface flux changes from   (Fig.6b) and  − (Fig.6f) are -3.9W/m 2 and -4.6W/m 2 respectively.Since the total energy tendency is almost zero over this area, the corresponding changes in vertical flux divergence (equal to net surface flux minus net TOA flux; Fig.6c) over this area are -1.8W/m 2 and -2.5W/m 2 respectively.The negative signs indicate that vertical flux divergence decreased and consequently divergence of horizontal energy transports increased in the 2001-2008 period compared to the 1986-2000 mean (compare equation 6), so both changes in TOA fluxes and atmospheric energy transport contribute roughly equally to the reduced downward surface fluxes over the eastern tropical Pacific from these two mass adjusted data sets.For   (Fig. 6d) the mean change in surface flux over this area is about -0.5W/m 2 and the corresponding mean change in vertical flux divergence (Fig. 6e) is about 1.6W/m 2 which is opposite to the mean changes in vertical flux divergence of   and  − , implying that increased horizontal energy transport into the east Pacific region offsets much of the reduction in TOA downward fluxes leading to a smaller change in surface fluxes in this case.The net surface flux change obtained from   is weaker than those obtained from   and  − , since   and  − are computed from analysed state quantities they are considered more realistic than   which is computed from model forecasts.Changes in TOA fluxes are about -0.5W/m 2 for UPSCALE and AMIP5 data (Fig. 6h,k).The changes at the surface (Fig. 6i,l) are 2.2W/m 2 and 3.3W/m 2 and the corresponding mean divergence changes of horizontal energy transport (Fig. 6j,m) are 2.7W/m 2 and 3.8W/m 2 , respectively, implying that increased horizontal energy transport by the atmosphere into the region dominate the simulated changes in the surface fluxes.The divergence difference over the eastern tropical Pacific between the mass adjusted data and those from model simulations requires further study.
et al.[2015]  highlighted the decline of TOA net downward radiation flux from 1999-2005 which potentially contributed to the recent warming slowdown.Consistent withSmith et al. [2015], similar calculations of two five year means centred at 1999 and 2005 from net downward surface energy fluxes show declines of 0.31 Wm -2 (reconstruction), 0.51 Wm -2 (UPSCALE), 0.07 Wm -2 (AMIP5) and 0.26 Wm -2 (ERA20C).The differences between flux changes at TOA and surface (Fig.6h-k) include the total energy tendency and divergence.The patterns are very similar to those surface changes, implies the atmospheric energy divergence is the dominant factor affecting the surface flux changes, since both changes of TOA flux and atmospheric energy tendency are relatively small.The changes of northward total meridional energy transport calculated from   ,   and  − are also plotted in Fig.4b.Energy transports from mass corrected divergences show the increase of northward transport in the northern hemisphere, but the energy transport from   shows a decrease.It is mixed in the south hemisphere where transport derived from   displays a small energy transport while both calculation from   and  − indicate an increase of poleward energy transport between 10─55 o S and 15─70 o S. The effect of the temporal

Fig. 1 .
Fig. 1.Schematic illustrating the energy flow terms used in the estimation of surface energy flux over land and ocean.

Fig. 2 .
Fig. 2. Deseasonalised anomaly (relative to the 2001-2005 period) time series of mean net downward radiation fluxes at TOA over (a) the globe, (b) the global ocean and (c) the global land, for data sets of AMIP5, ERAINT, WFOV,   and UPSCALE.Shaded areas of AMIP5 are sixteen member ensemble mean ± 1 standard deviation.All lines are 6 month running mean.The WFOV anomaly (60 o S-60 o N) is relative to 1985-1999 period, its line is three data points (three 72 day means) running mean and is adjusted to match   .The y-axis unit is W/m 2 on the left and PW on the right.(d) is the multi-annual (2001-2005) mean from   .The area mean (W/m 2 ) is displayed in the zonal mean plot.

Fig. 3 .
Fig. 3. (a) Multi-annual (2001-2005) mean net downward energy fluxes (in W/m 2 ) at surface from   .Zonal mean variations from AMIP5,   ,   , ERAINT, ERA20C, UPSCALE and  − are in the lower panel for (b) the globe, (c) the global ocean and (d) the global land, respectively.Shaded areas of AMIP5 are sixteen member ensemble mean ±1 standard deviation.The area mean is displayed in the zonal mean plot.

Fig. 6 .Fig 7 .
Fig. 6.Change in net energy fluxes (W/m 2 ,2001-2008 minus 1986-2000)  at TOA (left column), at surface (middle column) and the difference (right column) between fluxes at surface and TOA from reconstructions (  ,   and  − ), UPSCALE and AMIP5 data sets.a-c show the reconstruction based onAllan et al. [2014] at the TOA and the mass correction method using ERA Interim data, d-e are based on the residual method using ERA Interim data, f-g show the estimates from the mass correction method using MERRA reanalysis data, h-j are from the 5 ensemble mean of the UPSCALE simulations and k-m are the 16 ensemble member mean from the AMIP simulations.The marked area in b,d and f is from 20 o N-20 o S and 210 o E to the west coast of central America.
The total energy input to the atmosphere   =   −   where   is the net downward radiation flux (difference between the absorbed solar radiation and the outgoing longwave radiation) at TOA and

Table 1 .
Data sets and their properties o  Mizielinski et al. [2014]