# An assessment of air-sea heat fluxes from ocean and coupled reanalyses

Maria Valdivieso, Keith Haines, Magdalena Balmaseda, You-Soon Chang, Marie Drevillon, Nicolas Ferry, Yosuke Fujii, Armin Köhl, Andrea Storto, Takahiro Toyoda, et al.

► To cite this version:

Maria Valdivieso, Keith Haines, Magdalena Balmaseda, You-Soon Chang, Marie Drevillon, et al. An assessment of air-sea heat fluxes from ocean and coupled reanalyses. Climate Dynamics, 2017, 49 (3), pp.983-1008. 10.1007/s00382-015-2843-3. hal-03003585

HAL Id: hal-03003585

https://hal.science/hal-03003585v1

Submitted on 13 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# An assessment of air–sea heat fluxes from ocean and coupled reanalyses

Maria Valdivieso<fnref n="1" /> · Keith Haines<fnref n="1" /> · Magdalena Balmaseda<fnref n="2" /> · You-Soon Chang<fnref n="3" />,<fnref n="4" /> · Marie Drevillon<fnref n="5" /> · Nicolas Ferry<fnref n="5" /> · Yosuke Fujii<fnref n="6" /> · Armin Köhl<fnref n="7" /> · Andrea Storto<fnref n="8" /> · Takahiro Toyoda<fnref n="6" /> · Xiaochun Wang<fnref n="9" /> · Jennifer Waters<fnref n="10" /> · Yan Xue<fnref n="11" /> · Yonghong Yin<fnref n="12" /> · Bernard Barnier<fnref n="13" /> · Fabrice Hernandez<fnref n="5" />,<fnref n="14" /> · Arun Kumar<fnref n="11" /> · Tong Lee<fnref n="15" /> · Simona Masina<fnref n="8" /> · K. Andrew Peterson<fnref n="10" />

Received: 10 February 2015 / Accepted: 15 September 2015 / Published online: 5 October 2015
© The Author(s) 2015. This article is published with open access at Springerlink.com

**Abstract** Sixteen monthly air–sea heat flux products from global ocean/coupled reanalyses are compared over 1993–2009 as part of the Ocean Reanalysis Intercomparison Project (ORA-IP). Objectives include assessing the global heat closure, the consistency of temporal variability, comparison with other flux products, and documenting errors against in situ flux measurements at a number of OceanSITES moorings. The ensemble of 16 ORA-IP flux estimates has a global positive bias over 1993–2009 of 4.2 ± 1.1 W m⁻<fnref n="2" />. Residual heat gain (i.e., surface flux + assimilation increments) is reduced to a small positive imbalance (typically, +1–2 W m⁻<fnref n="2" />). This compensation between surface fluxes

and assimilation increments is concentrated in the upper 100 m. Implied steady meridional heat transports also improve by including assimilation sources, except near the equator. The ensemble spread in surface heat fluxes is dominated by turbulent fluxes (>40 W m⁻<fnref n="2" /> over the western boundary currents). The mean seasonal cycle is highly consistent, with variability between products mostly <10 W m⁻<fnref n="2" />. The interannual variability has consistent signal-to-noise ratio (~2) throughout the equatorial Pacific, reflecting ENSO variability. Comparisons at tropical buoy sites (10°S–15°N) over 2007–2009 showed too little ocean heat gain (i.e., flux into the ocean) in ORA-IP (up to 1/3 smaller than buoy measurements) primarily due to latent heat flux errors in ORA-IP. Comparisons with the Stratus buoy (20°S, 85°W) over a longer period, 2001–2009, also show the ORA-IP ensemble has 16 W m⁻<fnref n="2" /> smaller net heat

This paper is a contribution to the special issue on ocean estimation from an ensemble of global ocean reanalyses, consisting of papers from the Ocean Reanalyses Intercomparison Project (ORA-IP), coordinated by CLIVAR-GSOP and GODAE OceanView. The special issue also contains specific studies using single reanalysis systems.

<figure><img src="image_3.png" /></figure>

<footnote>Maria Valdivieso
m.valdiviesodacosta@reading.ac.uk</footnote>

<footnote n="1">¹ National Centre for Earth Observation (NCEO), Department of Meteorology, University of Reading, Harry Pitt Building, Whiteknights, PO Box 238, Reading RG6 6AL, UK</footnote>

<footnote n="2">² European Centre for Medium-Range Weather Forecasts (ECMWF), Reading, UK</footnote>

<footnote n="3">³ Geophysical Fluid Dynamic Laboratory, National Oceanic and Atmospheric Administration (GFDL/NOAA), Princeton, NJ, USA</footnote>

<footnote n="4">⁴ Department of Earth Science, Kongju National University, Kongju, South Korea</footnote>

<footnote n="5">⁵ Mercator-Océan, Toulouse, France</footnote>

<footnote n="6">⁶ Meteorological Research Institute (MRI), Japan Meteorological Agency (JMA), Tsukuba, Japan</footnote>

<footnote n="7">⁷ University of Hamburg, Hamburg, Germany</footnote>

<footnote n="8">⁸ Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC), Bologna, Italy</footnote>

<footnote n="9">⁹ Joint Institute for Regional Earth System Science and Engineering, University of California, Los Angeles, CA, USA</footnote>

<footnote n="10">¹⁰ Met Office, Exeter, UK</footnote>

<footnote n="11">¹¹ Climate Prediction Center, NOAA/NWS/NCEP, Camp Springs, MD, USA</footnote>

<footnote n="12">¹² Centre for Australia Weather and Climate Research, Bureau of Meteorology (BOM), Melbourne, Australia</footnote>

<footnote n="13">¹³ Laboratoire de Glaciologie et Géophysique de l'Environnement (LGGE), Centre National de Recherche Scientifique (CNRS), Grenoble, France</footnote>

<footnote n="14">¹⁴ Institut de Recherche pour le Developpement (IRD), Toulouse, France</footnote>

<footnote n="15">¹⁵ Jet Propulsion Laboratory (JPL), California Institute of Technology, Pasadena, CA, USA</footnote>
gain, nearly all of which is due to too much latent cooling caused by differences in surface winds imposed in ORA-IP.

**Keywords** Surface heat fluxes · Assimilation fluxes · Flux variability · Flux comparisons with in situ buoy flux data · Ocean and coupled reanalyses

# 1 Introduction

Surface heat fluxes over the ocean form a key component of the Earth's energy budget (Peixoto and Oort 1992) as the oceans comprise the main heat storage reservoir in the climate system and so determining how and where the oceans are warming up should provide important constraints for IPCC-class coupled climate models (e.g., Palmer and McNeall 2014). The air-sea heat fluxes are also needed to provide atmospheric forcing fields for ocean-only models (e.g., the Coordinated Ocean-ice Reference Experiments (COREs) documented in Griffies et al. 2009) and to assess heat budgets and the implied meridional transports of heat in comparison with direct oceanic transport estimates based on hydrographic section data (e.g., Bryden and Imawaki 2001; Macdonald and Baringer 2013). There has been increased interest recently in trying to improve air-sea heat flux estimates because of the availability of new satellite radiation data at the top of the atmosphere, e.g., from Clouds and Earth's Radiant Energy Systems (CERES; Loeb et al. 2009), and the availability of the Argo float network since the early 2000s for measuring changes in the upper 2000 m ocean heat content. A review of the current state of the art and the potential goals for the future can be found in WCRP (2012), Josey et al. (2013), Yu et al. (2013) and von Schuckmann et al. (2015).

Air sea fluxes are notoriously difficult to determine on large space (basin-scale) and timescales (interannual-to-decadal) mainly due to sensitivity of the fluxes to parameterizations of boundary layer processes, cloud radiative feedbacks and to the highly variable wind speed and sea state conditions; see WGASF (2000). Global products based on locally estimated quantities such as the National Oceanography Centre (NOC) flux products (Josey et al. 1999; Berry and Kent 2009), which used ship observations from the International Comprehensive Ocean-Atmosphere Data Set (ICOADS), often show unrealistically large heat flux biases (+15–30 W m⁻²) when integrated over global scales. Some atmospheric reanalysis products, such as the NCEP/NCAR reanalysis R1 (Kальнай et al. 1996; referred to as NCEP-R1) and the NCEP/DOE reanalysis R2 (Kanamitsu et al. 2002; referred to as NCEP-R2), have global budgets that are close to being balanced (~3 W m⁻²). However, others such as the ECMWF ERA-Interim reanalysis (Dee et al. 2011; referred to as ERAi), or the NASA

MERRA reanalysis (Rienecker et al. 2011), still have significant global heat budget imbalances (+11 W m⁻² for ERAi and +21 W m⁻² for MERRA, see Josey et al. 2013, Table 5.1). Furthermore, they cannot reliably reproduce surface fluxes that directly depend on clouds (such as radiation and precipitation) because cloud observations are not assimilated, and the turbulent fluxes also suffer from the absence of coupled feedbacks with the ocean due to the fixed surface temperature boundary conditions that do not allow ocean temperatures to respond, e.g., underneath midlatitude or tropical storms.

Josey and Smith (2006) argued that progress in flux field development requires careful evaluation of global flux products against independent observations, particularly from surface flux buoys. Many flux products now exist including some based on satellite data, such as the Japanese Ocean Flux data sets with Use of Remote sensing Observations (J-OFURO; Kubota et al. 2002), the Hamburg Ocean Atmosphere Parameters and Fluxes from Satellite Data (HOAPS, Andersson et al. 2010) and the IFREMER turbulent flux product documented in Bentamy et al. (2013). Examples of hybrid products (a combination of atmospheric reanalysis and satellite measurements) adjusted with satellite and/or in situ measurements include the Common Ocean-ice Reference Experiments (CORE.2) dataset (Large and Yeager 2009) and the WHOI Objectively Analyzed Ocean-Atmosphere Fluxes (OAFlux) product (Yu et al. 2008), which is the product that has been most completely evaluated against in situ flux measurements so far.

This paper is a first review of surface heat fluxes from an ocean reanalysis perspective, using the joint CLIVAR-GSOP/GODAE OceanView Ocean Reanalysis Intercomparison Project (ORA-IP) datasets (Balmaseda et al. 2015 and CLIVAR Exchanges no. 64). These ocean reanalyses (hereafter the ORA-IP products) are examples of ocean data assimilated models that are actively being used either for climate monitoring studies, e.g., ocean heat content (Xue et al. 2012; Balmaseda et al. 2013b) or steric sea level (Storto et al. 2015) variability, or for operational ocean forecasting (Lellouche et al. 2013; Blockley et al. 2014). Air-sea heat fluxes are made up of short and longwave radiation terms, along with turbulent fluxes for heat (sensible and latent) computed from bulk formulae, with both the outgoing longwave radiation (computed using the Stephan-Boltzmann Law) and the turbulent fluxes depending sensitively on the sea surface temperature (SST). For ORA-IP, air-sea fluxes are then generally computed based on prescribed atmospheric states from either atmospheric reanalyses or a blend of atmospheric reanalysis and satellite observations. However, because the SSTs are influenced by the near surface data being assimilated and used to constrain the different ORA-IP products, the results from ORA-IP still develop a range of estimated air-sea fluxes.

<figure><img src="image_4.png" /></figure>
This comparison examines both net surface heat fluxes and the individual flux components, using the ORA-IP ensemble means and spreads for years 1993–2009. The analysis is limited both by the variety of forcing design across the different products (details provided in Sect. 2), and the availability of output data from the products, therefore it is very hard to account for all the differences. A broader discussion of other ways forward is presented in the Sect. 6 at the end of the paper.

The paper outline is as follows. A summary of the forcing methods used for the ocean reanalyses, along with details of the assimilation increment diagnostics and additional flux data sets used in the comparison, are given in Sect. 2. Section 3 looks at the global time-mean heat budgets and compares these with global flux products from a variety of sources, including ship- and satellite-based products and atmospheric reanalysis fluxes. It also looks at the implied meridional heat transports and the role of the assimilation increments in closing the mean ocean heat budget. Section 4 looks at the temporal variability of the surface heat fluxes at seasonal and interannual timescales based on an ORA-IP ensemble of surface heat flux estimates. Section 5 makes comparisons of monthly mean ORA-IP surface heat fluxes with in situ flux data taken from the operational tropical moored buoy arrays and one OceanSITES station (Stratus) in the eastern South Pacific, where longer flux records exist. Section 6 provides a summary and further discussion.

## 2 Heat flux products

### 2.1 Surface heat fluxes in ORA-IP

Sixteen monthly surface heat flux products originating from ocean and coupled reanalyses have been compared over the 17-year period (1993–2009). These reanalyses have each been run with different model configurations, data assimilation systems and observational data sources. Table 1 summarises the ocean/coupled model frameworks and the choices made by each of the reanalysis for computing surface boundary forcing, including the atmospheric data sets, bulk formulae and SST observational products. Details beyond those provided in this table can be found in the cited references.

There are three main approaches that have been used to generate air-sea fluxes in ORA-IP: (1) “Flux Correction”, (2) “Bulk Flux Forcing”, and (3) “Coupled Model Fluxes”, although the variety of treatments within these classes is great and later results cannot be easily distinguished according to these classifications.

* “Flux Correction” for PEODAS, ORAS4 and GODAS

With this approach, the surface fluxes (for momentum, heat and freshwater) from an atmospheric reanalysis product are applied directly to the ocean surface, along with a surface relaxation of SST towards an observational product to prevent model drift. In ORA-IP, this SST restoring is often applied with rather short damping timescales (typically, 1–5 days) as a way of assimilating SST satellite gridded data into the models (see seventh column of Table 1).

* “Bulk Flux Forcing” for MOVECORE, MOVEG2, GECCO2, ECCOv4, CGLORS05v3, UR025.3, UR025.4, GloSea5, GLORYS2v1 and GLORYS2v3

In this approach, the turbulent fluxes (for heat, water and momentum) are derived from bulk formula using a prescribed atmospheric state and the model’s SST, which may also be affected by data assimilation (see Table 1). Nine of the reanalyses employ the CORE bulk formulae described in Large and Yeager (2004, 2009) applied using atmospheric data from either an atmospheric reanalysis (ERAi, NCEP-R1 or JRA-55) or a blend of reanalysis data and satellite observations (CORE.2). The near surface atmospheric variables at high temporal resolution (6- or 3-h) are adjusted from 2 to 10 m, following the Monin–Obukhov similarity parameterisation (Large and Yeager 2004; Eqs. 9b-c), and then combined with the model’s SST (model top-level potential temperature) and surface currents, to compute the turbulent fluxes at each model time step. The chosen atmospheric state also includes precipitation and runoff data (not discussed here), and radiative (downward shortwave and longwave) fluxes.

Some products (MOVEG2, CGLORS05v3 and GLORYS2v1/v3) applied a priori adjustments to atmospheric reanalysis data such as radiation and precipitation, to prevent biases associated with cloud parametrizations (Kallberg 1987). For example, CGLORS05v3 and GLORYS2v1/v3 applied corrections to radiative (shortwave and longwave) heat fluxes from the ERAi product by means of a large-scale climatological correction coefficient derived by the Global Energy and Water Cycle Experiment (GEWEX) project, following the method described by Dussin and Barnier (2013). Furthermore, two products use tuned Bulk Flux Forcing from the adjoint method (GECCO2, ECCOv4), where the surface fluxes are part of the control vector of the optimization problem (see Köhl 2015 for further details). GECCO2 fluxes are computed from Large and Pond (1981) bulk formula applied to NCEP-R1 fields, whereas ECCOv4 uses CORE bulk formula with all a priori forcing data taken from the ERAi product.

* “Coupled Model Fluxes” for MOVE-C, CFSR and ECDA
<table><caption>**Table 1** Summary of the surface heat flux datasets from ocean or coupled reanalyses (ORA-IP products) used in the comparison study</caption><thead><tr><th>ORA-IP<br>Data sets</th><th>Model</th><th>Data assimilation</th><th>Forcing method</th><th>Atmospheric<br>datasets</th><th>Bulk formula</th><th>SST</th><th>Period</th><th>References</th></tr></thead><tbody><tr><td>BOM PEODAS</td><td>MOM2<br>2° × (0.5°–1.5°)</td><td>EnKF (T, S, SST)</td><td>Flux correction</td><td>ERA40 until 2001,<br>NCEP-R2 after-<br>wards</td><td>—</td><td>Relaxation to weekly<br>Olv2SST (1 day)</td><td>1980–2012</td><td>Yin et al. (2011)</td></tr><tr><td>ECMWF ORAS4</td><td>NEMO<br>1° × (0.3°–1°)</td><td>3DVAR (T, S, SLA)</td><td>Flux correction</td><td>ERA40 until 1988,<br>ERAi afterwards</td><td>—</td><td>Relaxation to weekly<br>Olv2SST (3 days)</td><td>1960–2009</td><td>Balmaseda et al.<br>(2013a)</td></tr><tr><td>MRI/JMA<br>MOVECORE</td><td>MRI.COM<br>1° × 0.5°</td><td>3DVAR (T, S, SLA)</td><td>Bulk flux forcing</td><td>Daily (Qsw, Qlw)<br>and 6-h (u,v,tair,q)<br>at 10 m from<br>CORE2</td><td>Large and Yeager<br>(2004)</td><td>SST neither assimilated<br>nor relaxed. Assimi-<br>lates T-profiles from<br>WOD09/GTSPP</td><td>1948–2007</td><td>Fujii et al. (2015)</td></tr><tr><td>MRI/JMA<br>MOVEG2</td><td>MRLCOM<br>1° × (0.3°–0.5°)</td><td>3DVAR<br>(T, S, SST, SLA)</td><td>Bulk flux forcing</td><td>6-h (Qsw, Qlw, 10 m<br>(u,v), 2 m (tair,q))<br>from JRA-55, inc.<br>Qsw, Qlw correc-<br>tions</td><td>Large and Yeager<br>(2004)</td><td>Assimilates daily<br>MGDSST (operational<br>SST analysis by JMA)</td><td>1993–2012</td><td>Toyoda et al. (2013)</td></tr><tr><td>U. Hamburg<br>GECCO2</td><td>MIT<br>1° × (0.3°–1°)</td><td>4DVAR<br>(T, S, SST, SLA)</td><td>Bulk flux forcing</td><td>6-h NCEP-R1 fields<br>adjusted by the<br>adjoint method</td><td>Large and Pond<br>(1981), with<br>transfer coef-<br>ficients after<br>Large and<br>Yeager (2004)</td><td>Assimilates<br>AMSR/E + HadISST</td><td>1948–2011</td><td>Köhl (2015)</td></tr><tr><td>ECCOv4</td><td>MIT<br>1° × (0.3°–1°)</td><td>4DVAR<br>(T, S, SST, SLA)</td><td>Bulk flux forcing</td><td>6-h ERAi fields<br>adjusted by the<br>adjoint method</td><td>Large and Yeager<br>(2004)</td><td>Assimilates TMI-<br>AMSR/E + Olv2SST</td><td>1993–2010</td><td>Wunsch and Heim-<br>bach (2013)</td></tr><tr><td>NCEP GODAS</td><td>MOM3<br>1° × (0.3°–1°)</td><td>3DVAR (T, SLA)</td><td>Flux correction</td><td>6-h NCEP-R2</td><td>—</td><td>Relaxation to weekly<br>Olv2SST (5 days)</td><td>1980–2011</td><td>Behringer (2007)</td></tr><tr><td>CMCC<br>CGLORS05V3</td><td>NEMOv3.2<br>1/2° × 1/2°</td><td>3DVAR<br>(T, S, SST, SLA)</td><td>Bulk flux forcing</td><td>Daily (Qsw, Qlw)<br>and 3-h 10 m (u,v)<br>and 2 m (tair,q)<br>from ERAi, inc.<br>Qsw, Qlw correc-<br>tions</td><td>Large and Yeager<br>(2004)</td><td>Relaxation to daily<br>1/4° NOAA SST<br>(AVHRR + AMSR-E)<br>(12 days)</td><td>1990–2011</td><td>Storto et al. (2014)</td></tr><tr><td>U. Reading<br>UR025.3</td><td>NEMO v2.3<br>1/4° × 1/4°</td><td>OI (T, S)</td><td>Bulk flux forcing</td><td>Daily (Qsw, Qlw)<br>and 6-h 10 m (u,v)<br>and 2 m (tair,q)<br>from ERAi</td><td>Large and Yeager<br>(2004)</td><td>SST neither assimilated<br>nor relaxed. Assimi-<br>lates T-profiles from<br>EN3_v2a</td><td>1989–2010</td><td>Haines et al. (2012)</td></tr><tr><td>U. Reading<br>UR025.4</td><td>NEMO v3.2<br>1/4° × 1/4°</td><td>OI (T, S, SST, SLA,<br>IC)</td><td>Bulk flux forcing</td><td>Daily (Qsw, Qlw)<br>and 6-h 10 m (u,v)<br>and 2 m (tair,q)<br>from ERAi</td><td>Large and Yeager<br>(2004)</td><td>Assimilates satellite<br>SST (AVHRR + (A)<br>ATSR (1993–<br>2007) + GHRSSST<br>(2008–10)) + in situ<br>SST from ICOADS</td><td>1989–2010</td><td>Valdivieso et al.<br>(2014a, b)</td></tr></tbody></table>

<figure><img src="image_5.png" /></figure>
<table><caption>Table 1 continued</caption><thead><tr><th>ORA-IP<br>Data sets</th><th>Model</th><th>Data assimilation</th><th>Forcing method</th><th>Atmospheric<br>datasets</th><th>Bulk formula</th><th>SST</th><th>Period</th><th>References</th></tr></thead><tbody><tr><td><strong>Met Office- Glo-<br>Sea5</strong></td><td>NEMO v3.2<br>1/4° × 1/4°</td><td>3DVAR (T, S, SST,<br>SLA)</td><td>Bulk flux forcing</td><td>3-h (Qsw, Qlw, 10 m<br>(u,v), 2 m (tair,q))<br>from ERAi,</td><td>Large and Yeager<br>(2004)</td><td>Assimilates satellite<br>SST (AVHRR + (A)<br>ATSR (1993-<br>2007) + GHRSSST<br>(2008-10)) + in situ<br>SST from ICOADS</td><td>1993-2010</td><td>Blockley et al. (2014)</td></tr><tr><td><strong>Mercator GLO-<br>RYS2v1</strong></td><td>NEMO<br>1/4° × 1/4°</td><td>KF (T, S, SST, SLA)</td><td>Bulk flux forcing</td><td>3-h (Qsw, Qlw, 10 m<br>(u,v), 2 m (tair,q))<br>from ERAi, inc.<br>Qsw, Qlw correc-<br>tions</td><td>Large and Yeager<br>(2004)</td><td>Assimilates Reynolds<br>1/4° AVHRR-only SST</td><td>1993-2009</td><td>Ferry et al. (2012)</td></tr><tr><td><strong>Mercator GLO-<br>RYS2v3</strong></td><td>NEMO<br>1/4° × 1/4°</td><td>KF (T, S, SST, SLA,<br>IC)</td><td>Bulk flux forcing</td><td>3-h (Qsw, Qlw, 10 m<br>(u,v), 2 m (tair,q))<br>from ERAi, inc.<br>Qsw, Qlw correc-<br>tions</td><td>Large and Yeager<br>(2004)</td><td>Assimilates Reynolds<br>1/4° AVHRR-only SST</td><td>1993-2011</td><td>Ferry et al. (2012)</td></tr><tr><td><strong>MRI/JMA<br>MOVE-C</strong></td><td>JMA/MRI.COM<br>1° × (0.3° -1°)</td><td>3DVAR (T, S, SST,<br>SLA)</td><td>Coupled model fluxes</td><td>JMA</td><td>—</td><td>Assimilates daily COBE-<br>SST data compiled in<br>JMA</td><td>1993-2011</td><td>Fujii et al. (2009)</td></tr><tr><td><strong>NCEP CFSR</strong></td><td>CFSRv2/MOM4<br>0.5° x (0.25° -0.5°)</td><td>3DVAR (T)</td><td>Coupled model fluxes</td><td>CFSR</td><td>—</td><td>Assimilates Had-<br>ISST + daily<br>1/4° NOAA SST<br>(AVHRR + AMSR-E)</td><td>1980-2011</td><td>Xue et al. (2011)</td></tr><tr><td><strong>GFDL ECDA</strong></td><td>CM2.1/MOM4<br>1° × (0.3°-1°)</td><td>EnKF (T, S, SST)</td><td>Coupled model fluxes</td><td>CM2.1 with 6-h<br>tair, u and v from<br>NCEP-R2 assimi-<br>lated in ECDA</td><td>—</td><td>Assimilates Had-<br>ISST + daily<br>1/4° NOAA SST<br>(AVHRR + AMSR-E)</td><td>1993-2011</td><td>Chang et al. (2013)</td></tr></tbody></table>

Products listed in bold have individual heat flux components as well. All ORA-IP datasets are available on a common WOA 1° by 1° grid on a monthly mean basis. Abbreviations in the fifth column are: tair: air temperature; q: air humidity; u, zonal wind component; v, meridional wind component; Qsw, shortwave radiative flux; Qlw, longwave radiative flux
<table><caption>**Table 2** Additional surface heat flux data used in this comparison study</caption><thead><tr><th>Type</th><th>Data set</th><th>Resolution</th><th>Period</th><th>References</th></tr></thead><tbody><tr><td>Ship-based</td><td>NOC2.0</td><td>Monthly, 1°</td><td>1973–2009</td><td>Berry and Kent (2009)</td></tr><tr><td rowspan="4">Satellite-based</td><td>CERES</td><td>Monthly, 1°</td><td>2000–</td><td>Loeb et al. (2009)</td></tr><tr><td>ISCCP-FD</td><td>3 h, 2.5°</td><td>1984–2009</td><td>Zhang et al. (2004)</td></tr><tr><td>J-OFURO</td><td>Daily, 1°</td><td>1988–2008</td><td>Kubota et al. (2002)</td></tr><tr><td>HOAPS3.2</td><td>Monthly, 0.5°</td><td>1987–2008</td><td>Andersson et al. (2010)</td></tr><tr><td rowspan="4">Atmospheric reanalysis</td><td>ERA-Interim</td><td>6 h, T255</td><td>1979–</td><td>Dee et al. (2011)</td></tr><tr><td>JRA-55</td><td>Daily, 1.25°</td><td>1958–</td><td>Kobayashi et al. (2015)</td></tr><tr><td>MERRA</td><td>Hourly, 0.5°</td><td>1979–</td><td>Rienecker et al. (2011)</td></tr><tr><td>NCEP-R2</td><td>Hourly, T62</td><td>1979–</td><td>Kanamitsu et al. (2002)</td></tr><tr><td rowspan="3">Hybrid</td><td>CORE.2</td><td>Monthly, 1°</td><td>1948–2006</td><td>Large and Yeager (2009)</td></tr><tr><td>TOA CERES/ERAi Div.</td><td>Monthly, 1°</td><td>1984–</td><td>Liu et al. (2015)</td></tr><tr><td>OAFlux</td><td>Daily, 1°</td><td>1983–</td><td>Yu et al. (2008)</td></tr><tr><td rowspan="4">Buoy</td><td>TAO/TRITON</td><td>Monthly, Tropical Pacific</td><td>2007–2009</td><td>McPhaden et al. (1998)</td></tr><tr><td>RAMA</td><td>Monthly, 15°N90°E</td><td>2007</td><td>McPhaden et al. (2009)</td></tr><tr><td>PIRATA</td><td>Monthly, Tropical Atlantic</td><td>2007–2009</td><td>Servain et al. (1998), Bourlès et al. (2008)</td></tr><tr><td>WHOI Stratus</td><td>Monthly, 19.9°S, 85.3°W</td><td>2001–2009</td><td>Weller (2015)</td></tr></tbody></table>

These three products come from coupled models where the surface fluxes are a resolved part of the overall solution (ocean, atmosphere, land and sea ice states) to the coupled reanalysis (details provided in the cited references in Table 1). Although this approach should give more internally consistent flux treatments, these products may also be more weakly constrained by observational data.

## 2.2 Conservation of heat in ORA-IP

The majority of the ORA-IP products (except ECCOv4 and GECCO2) are affected by interior sources and sinks of heat associated with the temperature assimilation using sequential filtering schemes. This is also true for volume and salt conservation (see Schiller et al. 2013). When diagnosing heat budgets these additional sources/sinks of heat should be taken into account. To do this, we partition the total flux into the net heat flux, Qnet at the surface, plus a term arising from changes in temperature within the fluid column associated with assimilation increments for temperature, Qassim (referred to as the assimilation heat flux in the following):

$$Q_{tot} = Q_{net} + Q_{assim}, \qquad (1)$$

with $Q_{net} = Q_{sw} + Q_{lw} + Q_{lat} + Q_{sen}$ (the sum of the short-wave, longwave, latent and sensible heat fluxes), and Qassim determined (in heat flux units; W m⁻²) according to:

$$Q_{assim} = C_p \rho_0 \int_{-H}^{0} \frac{(T_{obs} - T_{model})}{\tau} dz, \qquad (2)$$

where the vertical integral of the local temperature increment ($T_{obs} - T_{model}$) extends over the full ocean column from the bottom $z = -H$ to the surface $z = 0$, and $\tau$ is the data assimilation time window (in sec), and $\rho_o = 1035$ kg m⁻³ and $C_p = 4000$ kg⁻¹ J K⁻¹ are representative sea water density and specific heat capacities, respectively.

## 2.3 Other surface flux products used in this comparison

The ORA-IP products are compared with other global air-sea heat flux data from in situ (ship) observations, satellite data, atmospheric reanalyses, or hybrid products (a combination of atmospheric reanalysis and remote sensing products), and locally, with buoy flux data measured at moorings (limited in both time and space). All additional flux products are listed in Table 2. These other global flux products are not references, but they do contribute to assessment of the uncertainty in the context of evaluating ORA-IP, while the point comparisons against local buoy data allow the seasonal cycle and annual mean fluxes from ocean reanalysis products to be evaluated and calibrated.

Except for the buoy data (whose details are provided in Sect. 5), all datasets listed in Tables 1 and 2 are available on a monthly mean basis and have been interpolated to a common 1° by 1° grid. All quantitative comparisons in the following sections are performed with monthly mean data on this common grid.
# 3 Global time-mean surface heat fluxes

## 3.1 Spatial maps of Qnet

Global maps of the 17 year (1993–2009) time-mean net surface heat flux, Qnet, for 16 individual ORA-IP products are shown in Fig. 1. While the overall patterns are similar there are many notable differences in magnitude. Overall, the coupled products (14–16) appear more positive than other products over large areas (this is shown in Fig. 2 where products are ordered by their global means). In particular, ECDA shows additional warming in the southern hemisphere while CFSR has the strongest heating rates in the western tropical Pacific and north Indian ocean. Among the eight products (1–8) forced with ERAi fields, CGLORS05v3 (3) shows visibly more cooling, especially in the tropical oceans. The cause may be related to an over-correction of the radiative fluxes from the ERAi product in this analysis. It is also seen that high latitude regions are not available for ECCOv4 (2) or GODAS (10) products, and in PEODAS (11) the high latitudes are clearly anomalous compared to other products.

## 3.2 Global heat closure

Figure 2 shows on the left global time-mean (1993–2009) total, net surface and assimilation heat fluxes from all ocean products, compared on the right to the global means of the other flux products listed in Table 2. The products (labelled on the x-axis of Fig. 2) have been ordered according to their global surface imbalance, from positive to negative. The 16-member ensemble mean of all ocean reanalysis fluxes is also shown (grey bar). A common global mask (all cells that are not flagged as land on the common grid) has been used for all estimates in Fig. 2 as this is the domain to which the whole Earth’s energy budget closure constraint applies (see e.g., Josey et al. 2013, Table 5.1).

For the ensemble of 16 ORA-IP products the 17-year average (±interannual STD) is $4.2 \pm 1.1 \text{ W m}^{-2}$, which increases to $6.8 \pm 1.3 \text{ W m}^{-2}$ when averaged ±60° latitude to avoid areas with some missing data (seen in Fig. 1; other models also have very crude treatments for sea-ice). The $2.6 \text{ W m}^{-2}$ difference in ocean heat gain reflects the high latitude cooling regions and equates to a poleward heat transport of 0.32 PW across 60°N into the Arctic, comparable to observational estimates of 0.28 PW across 55°N in the Atlantic (Bacon 1997). However, the $6.8 \text{ W m}^{-2}$ of ocean heat gain over ±60° latitude would clearly imply a much larger poleward heat transport, but is instead balanced by assimilation terms and heat storage gains (see below).

The ocean reanalysis net surface heat flux products (Qnet, blue bars in Fig. 2) have global imbalances in surface heating ranging from $+13 \text{ W m}^{-2}$ (for coupled products ECDA, CFSR) to $-10 \text{ W m}^{-2}$ (CGLORS05v3), generally showing a similar level of closure to the atmospheric reanalysis products, between $+11.5 \text{ W m}^{-2}$ (MERRA) and $-15 \text{ W m}^{-2}$ (JRA-55). The independent ship and satellite products (e.g., NOC2.0, ISCCP/J-OFURO and ISCCP/OAFlux) have larger imbalances, $+15-25 \text{ W m}^{-2}$, although these may be biased by the ISCCP radiation product.

The ocean reanalyses generally have a much smaller global heat budget residual when the data assimilation terms, Qassim defined in Eq. (2), are taken into account to produce an equivalent total heat source, $Qtot = \text{Qnet} + \text{Qassim}$ (green bars in Fig. 2). For seven of the products this total flux, Qtot, gives a heat gain between $0.5 \text{ and } 2 \text{ W m}^{-2}$, which is close to recent estimates of heat content change from combined XBT and Argo data, which give warming rates of $0.64 \pm 0.11 \text{ W m}^{-2}$ (0–700 m) between 1993 and 2008 (see Loeb et al. 2012 and Roemich et al. 2015). For most of the products the assimilation terms, Qassim (orange bars in Fig. 2), represent a global ocean heat loss, countering the heat gain through the resolved surface heat flux, and we will see below that this cancellation mostly takes place in the top 100 m.

## 3.3 Oceanic feedback on heat fluxes through SST

Figure 3a presents timeseries for the 60°S–60°N averaged sea-air temperature differences (upon which the turbulent heat fluxes sensitively depend) for six ORA-IP products forced by CORE Bulk Formula using ERAi surface fields. The ERAi line (bold dashed black) represents the difference between the air temperatures and the SSTs originally used as atmospheric lower boundary conditions (see Dee et al. 2011), and this is seen to be fairly steady ~1.1° C. All products show positive values reflecting the ocean being warmed by radiation and needing to cool through turbulent fluxes. Some results are relatively steady like ERAi, whereas others show significant discrepancies, for instance, UR025.4, GloSea5 and GLORYS2v1/v3 all show smaller values in earlier years suggesting that the assimilated satellite SSTs (a combination of AATSR, AVHRR and AMSRE data) are cooler than the Reynold’s SSTs used to determine ERAi air temperatures.

The net surface heat fluxes into the ocean (Qnet) from ERAi and from the six ORA-IP products, are shown in Fig. 3b. The SST rise seen in some of the products over years 1993–2009 (Fig. 3a) leads to a decrease in Qnet (i.e., an increase in ocean cooling dominated by latent heat), whereas the ERAi Qnet shows no obvious trend. The assimilation heat fluxes in the top 100 m for these
<figure><img src="image_6.png" /><figcaption>**Fig. 1** Global spatial distribution of net surface heat fluxes from ORA-IP products (listed in Table 1) averaged over their common years (1993–2009). MOVECORE (panel 12) was based on the CORE.2 fluxes which end in 2007, and so was extended through 2009 using the monthly CORE.2 climatology. Positive is heat flux into the ocean. Contours indicate the locations of the 17-year mean zero net heat fluxes. Colorbar units are in W m⁻²</figcaption></figure>

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 2** Global mean heat fluxes averaged over the 17-year period (1993–2009) along with their interannual STDs over this period. Sixteen ORA-IP Qnet (net surface heat flux) products (blue bars) along with the 16-member ensemble average (dark grey bar) are shown in comparison with other products derived from observations, atmospheric reanalyses or blended products (light grey bars) to the right-hand side, with the error bars representing interannual</figcaption></figure>

standard deviations. The products (labelled on the x-axis) have been ordered according to their global surface imbalance, from positive to negative. Nine ORA-IP products have assimilation heat fluxes, Qassim, defined in Eq. (2) (orange bars), along with total heat fluxes, Qtot = Qnet + Qassim (green bars). A common global ocean-land mask has been used for all estimates (see text for details). Positive is heat flux into the ocean. Units are in W m⁻²

six products (calculated from Eq. (2) with the vertical integral now extending from the surface down to 100 m) are shown in Fig. 3c. These are generally negative (cooling the upper 100 m of the ocean, consistent with the full depth increments from Fig. 2) and are also decreasing with time. For example, the surface heat flux change of more than 10 W m⁻² in GloSea5 is balanced by a decrease in the upper layer assimilation cooling by about the same amount. This is strong evidence that there is compensation between near surface assimilation increments (mostly resulting from SSTs) and the applied surface heat fluxes. This does suggest that assimilation increments may well contain useful information about surface heat flux biases.

## 3.4 Ensemble comparisons of mean surface heat fluxes

Figure 4 shows the ORA-IP ensemble mean and spread (STD) of the 17-year (1993–2009) time-mean net surface heat flux, Qnet, and also the spread of the net radiative Qrad (shortwave plus longwave), and turbulent Qtur (latent plus sensible) heat fluxes using all products in Table 1.

The ensemble mean Qnet (Fig. 4a) is in broad agreement with the climatological fluxes based on bulk formula applied to observations (e.g., Berry and Kent 2009), with the tropical oceans showing heat gain, and the subtropics and high latitudes heat loss, especially in the vicinity of the western boundary currents (WBCs). In most
<figure><img src="image_9.png" /><figcaption>**Fig. 3** Year-to-year variability of global (60°S–60°N) averaged: **a** Sea-air temperature (SST–Tair) from a subset of ORA-IP products forced by CORE Bulk Formula using ERAi surface fields in comparison with the ERAi product itself; **b** Net surface heat fluxes (the sum of shortwave, longwave, latent and sensible heat fluxes) from</figcaption></figure>

the ERAi product and as calculated from the six ORA-IP products; **c** Assimilation heat fluxes arising from changes in temperature associated with assimilation increments for temperature over the upper 100 m (calculated from Eq. (2)) with the vertical integral extending from the surface to 100 m)

areas the ensemble spread in Qnet (Fig. 4b) is dominated by turbulent heat fluxes (Fig. 4d), with the largest spreads (>40 W m⁻²) occurring over the northern hemisphere WBCs, as well as in the Southern Indian ocean south of 20°S. There is also significant spread contributed by net radiation (~25 W m⁻²), particularly around the cooler upwelling regions off the west coasts of continents

(Fig. 4c). If short and longwave radiation spreads are independently assessed they are found to partially compensate, being anti-correlated in most regions (not shown). These differences in shortwave radiative fluxes are related with problems simulating clouds in the atmospheric reanalysis products (ERAi, NCEP-R2 and JRA-55) and also to errors in the coupled models (e.g., the shortwave flux in CFSR is

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>

**Fig. 4** **a** Ensemble mean of the 17-year (1993–2009) time-mean net surface heat flux (Qnet) from the 16 ORA-IP products listed in Table 1. Ensemble spread (STD) of time-mean, **b** Qnet (the sum of the shortwave, longwave, latent and sensible heat fluxes), **c** net radiative Qrad (shortwave plus longwave) fluxes, **d** net turbulent Qtur (latent plus sensible) heat fluxes, and **e** the sum of spread in net radiative and turbulent heat fluxes (STD Qrad + STD Qtur). **f** Correlation (CORR) between ensemble member Qrad and Qtur fluxes over

<figure><img src="image_12.png" /></figure>

the same period. Positive values in Plate **a** indicate heat flux into the ocean (W m⁻²). The ensemble spread in Qrad (Plate **c**) is computed using 10 different shortwave radiation products (ERA40, ERAi, CORE.2, JRA-55, NCEP-R1, NCEP-R2, ERAi with corrections, MOVE-C, CFSR, and ECDA) and 16 different longwave radiation products (because these depend on the model’s SST). *Contours* of 15 W m⁻² (*black solid lines* on Plates **b–d**) indicate the global average of spread in long-term mean Qnet

much too strong in the western tropical Pacific compared to satellite radiation from ISCCP and CERES products). Clearly STD Qnet (Fig. 4b) ≪ [STD Qrad + STD Qtur] (Fig. 4e) indicating anti-correlation between Qrad and Qtur between members of the ensemble (Fig. 4f); standard deviations are shown rather than variances because the scales are easier to interpret. The exceptions being the area south of Japan where radiation dominates the spread in fluxes, perhaps due to difference in representing aerosols downstream of China, and in the southwest Indian ocean, where correlations between Qrad and Qtur are also positive.

## 3.5 Implied ocean heat transports

Zonal averages in the net surface heat fluxes can be used to infer a meridional ocean heat transport that is consistent with a steady state ocean circulation under these fluxes. Figure 5a, **b** show inferred meridional heat transports for the ensemble mean, and each of the ocean reanalyses. Figure 5a uses the resolved surface heat fluxes, i.e., $Qnet = Qsw + Qlw + Qlat + Qsen$, and Fig. 5b uses the total heat fluxes including the assimilation terms, i.e., $Qtot = Qnet + Qassim$ (Qassim defined in Eq. 2). The
<figure><img src="image_13.png" /><figcaption>**Fig. 5** a Global Meridional Heat Transport (MHT) inferred from integrated ORA-IP surface heat fluxes (starting from the south) and a steady state assumption (16 colour lines with the shading indicating ORA-IP ensemble spread). **b** MHT inferred from integrated surface heat fluxes adjusted by assimilation increments (nine products). The</figcaption></figure>

red symbols represent WOCE-based inverse model estimates at control sections from Ganachaud and Wunsch (2003) and Lumpkin and Speer (2007). Assimilation increments generally improve agreement with external transport estimates outside the tropical region ±10°. Positive numbers indicate northward transport. Units are in PW

fluxes for the reanalyses UR025.3/4 have been previously compared to the advective heat and freshwater transports in Haines et al. (2012) and Valdivieso et al. (2014a); however, the full reanalysis transports have not yet been compared for the full range of reanalyses shown here. In Fig. 5a, **b** all surface fluxes are integrated starting in the south (i.e., the Antarctic continent) and working northwards. At a number of latitudes an independent meridional heat transport is provided based on inverse modelling using hydrographic section data from Ganachaud and Wunsch (2003) and Lumpkin and Speer (2007).

Figure 5a shows a rapid divergence of results moving northwards through the southern ocean, with the largest global heat flux imbalance in ORA-IP, ~13 W m⁻² for CFSR and ECDA (Fig. 2), corresponding to an integrated transport discrepancy in the north ~+5 PW. The ensemble spread

from the ORA-IP products is shaded around the mean, and we see that this area grows rapidly in the southern ocean, remains steady in the southern subtropics, then grows again crossing the tropics, but remains fairly stable through the northern subtropical and higher latitudes. This suggests that the largest uncertainty in net surface heat fluxes occurs in the southern oceans and in the tropics, which is broadly consistent with the Qnet spread map in Fig. 4b.

When the assimilation terms are added in Fig. 5b, the reanalysis products are brought much closer to a steady balance, with most products showing less than 0.5 PW residual transports at 80°N, reflecting a global net flux imbalance now less than 1.5 W m⁻². The shaded area in Fig. 5b now represents the 9-member ORA-IP spread which is greatly reduced compared to Fig. 5a (the same 9 ORA-IP members on Fig. 5a have almost the same mean and spread as shown

<figure><img src="image_14.png" /></figure>
for all 16). However, in some of the products there are now highly anomalous implied ocean transports appearing near the equator. The anomalous transports are symmetric with strong values away from the equator, but converging just to the north and south, without obvious influence at higher latitudes (the shaded spread increases at the equator but is reduced again further north). The problem is most likely related to difficulties in producing a realistic equatorial dynamical balance in the models due to the lack of a pressure gradient bias correction in some of the data assimilation systems (see e.g., Bell et al. 2004; Balmaseda et al. 2007, for further details).

# 4 Evaluation of flux variability

In this section we look at the ORA-IP ensemble consistency of the net surface heat flux temporal variability, both for the seasonal cycle and on interannual timescales. We also evaluate regional co-variability of surface heat fluxes and SSTs from ORA-IP and other products on interannual time scales. Note that the comparisons hereafter are limited to the surface heat fluxes (Qnet) because the assimilation fluxes (Qassim defined in Eq. 2) are not available for some of the products.

## 4.1 Surface heat flux seasonal cycles

The consistency of the seasonal cycles among all ORA-IP Qnet (net surface heat flux) products in Table 1 is presented in Fig. 6. Figure 6a shows the standard deviation of the monthly mean (1993–2009) ensemble average Qnet (note that the mean flux biases are removed from each product). This increases away from the tropics and is largest around the northern WBCs where warm water is always present off the continental shelves, and where large heat losses in winter can be sustained when cold air flows off the continents. Figure 6b looks at the typical variability amongst seasonal cycles from each model based on monthly deviations from the ensemble mean seasonal cycle. Figure 6c, d show the same but for only the northern winter months (December–February) and southern winter months (June–August), respectively. The 10 W m⁻² contour is highlighted showing that the monthly mean variability amongst all the products generally differ by less than this over large areas. The northern western boundary currents, especially in the winter months, show the largest differences (>25 W m⁻²); monsoon upwelling areas off east and west Africa and the Arabian Peninsula also show large variability (due to short-wave and latent heat flux variability) in northern summers.

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 6** a Seasonal STD of Qnet (net surface heat flux) estimates (16 products) computed from the monthly climatology over 1993–2009 (the mean Qnet over the 17-year base period has been removed from each product). Ensemble spread of the 17-year mean seasonal</figcaption></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>cycle, b annual average (January–December), c Northern winter season (December–February), and d Northern summer season (June–August). The 10 W m⁻² contour is highlighted</figcaption></figure>
## 4.2 Interannual heat flux signals

To determine the robustness of the interannual net surface heat flux signals in the ensemble of products, Fig. 7 shows the Signal to Noise ratio for the interannual variability over the period 1993–2009. The Signal (Fig. 7a) is the temporal standard deviation of the ensemble annual mean Qnet anomalies from the 1993–2009 average, and the Noise (Fig. 7b) is the product standard deviation around the ensemble mean averaged over the same period. Two versions of the Signal to Noise ratio are shown; Fig. 7e using all the products, except PEODAS, and Fig. 7f using only 12 products, i.e., without the coupled reanalyses. The temporal STD of ensemble annual mean turbulent (latent plus

sensible) heat flux (Qlat + Qsen) and SST anomalies from 15 products except PEODAS over the same period, are also shown in Fig. 7c, d respectively, for comparison with Fig. 7a.

The strongest interannual signals in Qnet (Fig. 7a) are dominated by variability in latent and sensible heat loss from the ocean surface (Fig. 7c) and occur primarily in the El Niño-Southern Oscillation (ENSO) region in the tropical Pacific, where Qnet anomalies are strongly co-located with SST anomalies (Fig. 7d). Figure 7e shows interannual Signal to Noise Ratio up to 2 throughout the equatorial Pacific, reflecting the detection of ENSO, with the areas of detectable signal spreading to 20° N/S in the western Pacific. In the North Pacific, there is a suggestion of consistent

<figure><img src="image_19.png" /><figcaption>**Fig. 7** Interannual signals over 1993–2009 of a net surface heat fluxes, c turbulent (latent plus sensible) heat fluxes and d SSTs—all estimated from yearly anomalies relative to the 17-year (1993–2009) mean using all ORA-IP products, except PEODAS. b Noise estimated as the product STD around the ensemble mean averaged over</figcaption></figure>

the same period. e Signal to Noise Ratio for the ORA-IP Ensemble Qnet (15 products) and f Signal to Noise Ratio using only 12 products, excluding the coupled reanalyses. The *solid contours* on Plates e–f indicate the location of Signal/Noise ratios of 1

<figure><img src="image_20.png" /></figure>
interannual signal near the Gulf of Alaska, where val-
ues reach 1.2–1.3, that may be associated with the Pacific
Decadal Oscillation (PDO). However, in other areas of
large flux variability, such as the northern WBCs and their
extensions, there is little detectable interannual signal due
to large noise (>20 W m⁻²) among the products (Fig. 7b).
The coupled products contribute to a larger ensemble noise
in the tropics (particularly the CFSR; Fig. 1) and when
removed, the Signal/Noise in Qnet increases substantially,
especially in the western Pacific warm pool (Fig. 7f).

To develop a more detailed comparison of these signals,
Fig. 8 shows a set of longitude-time (Hovmöller) plots of
the monthly latent heat flux (Qlat) anomalies (removing
the 1993–2009 mean seasonal cycles) for the ENSO region
of the tropical Pacific (5°N–5°S, 130°E–80°W) from 1993
to 2009. Eleven of the ORA-IP products are represented
(wherever latent fluxes are available separately) along with
the 11-member ensemble mean, and eight other latent heat
flux products, including atmospheric reanalyses, ship- and
satellite-based or hybrid products.

Differences in the strength and location of the heat flux
anomalies associated with the El Niño 1997/1998 and La
Niña 1999/2000 are clearly seen in Fig. 8. During the build
up of El Niño in 1997 there is strong warming of the east-
ern equatorial cold tongue that coincides with the largest
increase of latent heat loss (Qlat > 0) around 220–260°E in
all products, except NOC2.0 which is too noisy (Fig. 8m).
The La Niña phase from 1999 to 2000 is marked by the
opposite behaviour, i.e., large negative Qlat anomalies
(reduced ocean heat loss) can be seen in most products,
responding to cool SSTs, but unlike the 1997/1998 flux
anomalies, they tend to spread eastward as their ampli-
tudes decrease. Among the ORA-IP products, the coupled
CFSR and ECDA are clear outliers, with CFSR showing
large flux anomalies in the central Pacific prior to 1997
and ECDA with positive/negative anomalies in the western
Pacific (130°–180°E) from 1993 to 2001, that are absent in
all other products. Most of the ERAi-forced ORA-IP prod-
ucts have spurious (non physical) positive flux anomalies
coincident with the Tropical Ocean–Atmosphere (TAO)
mooring lines, particularly from 1993 to 1995, that are not
in the other products, including in the original ERAi prod-
uct itself (Fig. 8q). Josey et al. (2014) noted this anomaly
pattern may be related to assimilation of near-surface data
from the TAO array in the ERAi reanalysis. The satellite
and hybrid products show quite consistent Qlat anomaly
patterns (Fig. 8n–p), whereas NCEP-R2 (Fig. 8r) has very
weak El Niño 1997/1998 and large anomalies in the central
Pacific after 2001.

We now look at the consistency in reproducing surface
flux anomalies during a particular anomalous year. Figure 9

shows ensemble ORA-IP SST (Fig. 9a) and Qnet (Fig. 9b)
anomalies in 2008 (relative to the period 2001–2009) in the
Pacific sector. The SST anomaly pattern in 2008 is associ-
ated with the Pacific Decadal Oscillation (PDO) (the most
negative since 1971; Peterson and Baringer 2009), with
negative anomalies along the west coast of North America
from Alaska to the equator, and positive anomalies in the
area to the west extending to up 30°N/S. Corresponding
net surface heat flux anomalies (positive Qnet corresponds
to oceanic heat gain) can be seen in the North Pacific up
to 40°N, with anti-correlation, i.e., increased Qnet (due to
reduction in latent and sensible heat loss) over negative
SST anomalies and reduced Qnet over positive SST anom-
alies, indicating where surface heat fluxes act to damp SST
anomalies that have been generated directly by the atmos-
phere. For comparison with the ORA-IP ensemble output,
Fig. 9c, d show the 2008 Qnet anomalies from the ERAi
product and the CERES radiation combined with the OAF-
lux dataset, respectively. Both the ORA-IP ensemble and
the ERAi reanalysis are capable of producing a very simi-
lar pattern of variability to that obtained from the combined
satellite CERES radiation and OAFlux fields for much of
the Pacific basin. This result is encouraging for studies that
rely mostly on the use of atmospheric reanalysis-based
radiation with potential biases with simulated clouds.

We also consider surface heat flux anomalies in 2009 in
the North Atlantic (Fig. 10) between July and December,
which are associated with a persistent negative phase of
the North Atlantic Oscillation (NAO) (Cayan 1992), which
started in July 2009 (Arndt et al. 2010). The NAO tripole
pattern is clearly seen in the Qnet anomaly pattern from the
ORA-IP ensemble (15 products except PEODAS) and the
CERES/OAFlux combined product (Fig. 10a, b, respec-
tively). Differences between the ensemble of ORA-IP and
CERES/OAFlux toward the mid-latitude western bounda-
ries suggest incorrect positioning of the Gulf Stream and
the North Atlantic Current, which is too zonal in some of
the models. The ERAi atmospheric analysis itself gives the
surface heat flux product in Fig. 10c, although this is influ-
enced by the fixed daily SST used at the surface boundary.
An alternative way of using the atmospheric reanalysis data
is to combine CERES TOA (top-of-atmosphere) fluxes
with a correction based on vertically-integrated ERAi heat
flux divergences in the atmosphere, to derive a net heat flux
at the ocean-atmosphere interface. Such a product is docu-
mented in Liu et al. (2015) and the results can be seen in
Fig. 10d. This flux product now contains short-scale anom-
alies associated with atmospheric winds that are not seen
in the original ERAi net heat flux anomalies (Fig. 10c),
but is otherwise consistent with the CERES/OAFlux
combination.
<figure><img src="image_21.png" /><figcaption>**Fig. 8** Monthly latent heat flux (Qlat) anomalies relative to the 17-year (1993–2009) base period in the Tropical Pacific (130°–280°E, meridionally integrated between 5°S and 5°N). Eleven of the ORA-IP products are represented along with the 11-member ensemble average (Plates a–l) in comparison with eight other latent heat flux products, including atmospheric reanalyses, ship and satellite derived products or blended products. *Positive anomalies* indicate more than normal oceanic heat loss due to increased latent heat flux</figcaption></figure>

<figure><img src="image_22.png" /></figure>
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

**Fig. 9** a Pacific basin-scale SST anomalies in 2008 (relative to 2001–2009) from the ensemble of ORA-IP (15 products, except PEODAS) and the corresponding anomalies in net surface heat fluxes, Qnet, as derived, b from the ORA-IP ensemble, c from the ERAi reanalysis and d from CERES radiation combined with OAF-lux turbulent fluxes. The solid contours represent the location of the

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

zero SST/Qnet anomalies. Positive Qnet anomalies represent more heat than normal going into the ocean in areas of mean net heat gain (tropics and eastern Pacific) or less heat than normal being lost from the ocean in areas of mean net heat loss (e.g., WBCs and their extensions), see Fig. 4a for the climatological location of the zero net heat flux

# 5 Ensemble comparisons with buoy flux data

Local evaluation of gridded surface flux products e.g., from models or satellite data, against in situ flux measurements is really only possible at a few locations where conditions controlling surface fluxes have been monitored by meteorological buoys for periods of a year or more to cover the seasonal cycle. In this section, we present a comparison between monthly mean ORA-IP surface heat fluxes, and the corresponding fluxes derived from buoy measurements taken from the operational tropical moored buoy arrays and the Woods Hole Oceanographic Institution (WHOI) Stratus buoy which is located in the eastern South Pacific. Such comparisons are very important to distinguish between the many available gridded flux products, as noted by Yu et al. (2013).

## 5.1 Comparison with tropical buoys

The comparisons begin with tropical buoys from the TAO array in the tropical Pacific (McPhaden et al. 1998), the RAMA array in the tropical Indian (McPhaden et al. 2009)

and the Pilot Research moored Array (PIRATA) in the tropical Atlantic (Servain et al. 1998; Bourlès et al. 2008). These buoy flux data are available through the OceanSITES project (http://www.pmel.noaa.gov/tao/oceansites/flux/main.html) and include radiative (shortwave and long-wave) and turbulent (latent and sensible) fluxes computed using the COARE3.0b flux algorithm (Fairall et al. 2003). Here we use eight buoy deployments within the 10°S–15°N latitude band for the 24 month period Jan 2008–Dec 2009 (except Jan–Dec 2007 for the PIRATA site 15°N38°N and Jan–Dec 2009 for the RAMA site 15°N90°E), when all four components of the heat flux are available. We note that long-moored timeseries are available at some locations (e.g., the TAO sites on the equator), but they are not continuously distributed in time. A model-data comparison for the full ORA-IP time frame 1993–2009 would require dealing with any gaps in the buoy data and interannual variations in surface heat flux computations associated with ENSO and is beyond the scope of this manuscript.

Figure 11a shows annual mean net heat flux and the individual flux components at the eight buoys, along with monthly standard deviations about the annual means
**Fig. 10** Net surface heat flux (Qnet) anomalies between July and December 2009 (computed as departures from 2001 to 2009 monthly means) in the North Atlantic sector: **a** ORA-IP ensemble using 15 products, except PEODAS, **b** CERES net radiation combined with OAFlux turbulent fluxes, **c** ERAi product (the sum of the shortwave, longwave, latent and sensible heat fluxes) and **d** top of atmosphere (TOA) CERES net radiation combined with ERAi transport divergences. The solid contours represent the location of the zero Qnet anomalies. *Positive Qnet anomalies* represent more heat than normal going into the ocean in areas of mean net heat gain (low latitudes regions south of 30°N) or less heat than normal being lost from the ocean in areas of mean net heat loss (high latitude regions north of 30°N)

<figure><img src="image_27.png" /></figure>

representing mainly the seasonal cycle. There is net heat flux into the ocean at all locations, with a maximum heating rate ~175 W m⁻² at the TAO site 0°N110°W in the eastern Pacific cold tongue. The net heat fluxes off the equator at 15°N and 10°S are considerably less than the net heat flux at the equator, reflecting wind-induced evaporation, particularly at the southeast PIRATA site 10°S10°W. The longwave and sensible heat fluxes vary rather little between buoys, and also have small seasonal variability, so that shortwave and latent heat fluxes are the dominant heat flux terms. Wang and McPhaden (2001, Appendix B) estimated the accuracy range for annual averaged flux components from TAO buoys (at 110°W, 140°W, 170°W and 165°E) to be: (1, 7–12, 6, 12–19 W m⁻²) for (sensible, latent, longwave, shortwave) fluxes, respectively. The net flux errors (assuming they are uncorrelated) are then 17, 23, 18, and 16 W m⁻² at 165°E, 170°W, 140°W and

110°W, respectively, always dominated by the shortwave errors.

Figure 11b shows (Ensemble ORA-IP–Buoy) differences using buoy annual heat fluxes and ensemble fluxes from the 11 products with all component fluxes available (listed in bold in Table 1), interpolated to each buoys location and sampling time period, with the STD now indicating spread among the ORA-IP products about the ensemble mean differences. “All Buoys” indicates average differences across all eight buoys.

Differences (biases) for nearly all flux components are negative indicating overestimation of ocean heat loss or underestimation of ocean heat gain. Shortwave biases, typically –6 ± 6 W m⁻², suggest the reanalysis downward shortwave fluxes are slightly too weak on average. Longwave net losses are 5 ± 3 W m⁻² too strong compared to the buoys indicating that the reanalyses have either a too
**Fig. 11** *a* Observed annual mean net surface heat fluxes and their individual flux components (i.e., shortwave, longwave, latent and sensible heat fluxes) at eight buoy locations of the operational tropical moored buoy arrays (10°S–15°N) over 2007–2009. The *error bars* represent monthly standard deviations. *Positive values* indicate heat flux into the ocean (W m⁻²). *b* Mean flux differences from buoy, (Ensemble ORA-IP-Buoy), derived from 11 ORA-IP products (listed in bold in Table 1), interpolated to each buoy location and averaged over the same period, with the *error bars* representing the spread among the various ORA-IP products. “All Buoys” indicates mean fluxes/ differences from buoy averaged across all eight buoys (see text for details)

<figure><img src="image_28.png" /></figure>

strong upward component or too weak downward compo-
nent. Both radiation biases contribute to cooling the ocean
model surface. Sensible heat flux biases are smaller lead-
ing to ~3 ± 1.5 W m⁻² additional cooling in the reanalysis
products. However, the latent flux biases dominate, with
~15 ± 9 W m⁻² of greater surface cooling in the ORA-
IP ensemble. The resulting net heat flux into the ocean
in ORA-IP is 10–41 W m⁻² (averaging at 29 W m⁻²) too
weak across all buoys, compared to the buoy average net
flux of ~87 W m⁻². None of the flux component errors are
correlated with the observed buoy fluxes themselves (as
seen from the distributions in Fig. 11a, b). However, flux
component errors are correlated between ensemble mem-
bers so that the range in net flux errors shown in Fig. 11b
are typically smaller than would be expected based on

errors in the components. The dominant effects are anti-
correlated shortwave and longwave flux errors and anti-cor-
related shortwave and latent flux errors.

Figure 12 shows a more detailed comparison in the western Pacific warm pool (165°E) and the equatorial cold tongue (110°W), with individual members of the ORA-IP ensemble and other products in Table 2, including the ICOADS-based product, five satellite-based or hybrid products and four atmospheric reanalysis products.

Negative offsets in nearly all components in all prod-
ucts indicate the excess cooling in ORA-IP, with biases
~25–40 W m⁻² in net flux not uncommon, well outside
of the TAO error estimates of ±17 W m⁻² at 165°E and
±15.6 W m⁻² at 110°W, also dashed in Fig. 12. Among the
ORA-IP products, MOVEG2 and the coupled CFSR and
<figure><img src="image_29.png" /><figcaption>**Fig. 12** Mean heat flux component differences from two TAO locations of the equatorial Pacific: **a** the western Pacific warm pool (165°E) and **b** the equatorial cold tongue (110°W) for the period</figcaption></figure>

2007–2009. The dashed vertical lines indicate the standard error for seasonal variations in net heat fluxes (Net) from TAO data based on Wang and McPhaden (2001, Appendix B)

ECDA reanalyses show the greatest agreement at 165°E, within error bounds of TAO data, but with bias compensation between radiative and latent flux components. Satellite radiation (ISCCP, CERES) combined with OAFlux, J-OFURO or HOAPS3.2 turbulent fluxes are all within the TAO error estimates, while the net flux bias in the ship-based NOC2.0 product is much larger, especially in the cold tongue region where difference reaches ~−60 W m⁻². These NOC2.0 differences come from the turbulent (latent and sensible) components, perhaps suggesting different biases in the bulk variables in the flux computation. For the atmospheric reanalysis products, differences of 35–75 W m⁻² less heat into the ocean are dominated by latent and shortwave fluxes, similar to the ORA-IP biases. The NCEP-R2 product underestimates shortwave radiation by about 50 and 65 W m⁻² compared with TAO. As was pointed out by Wang and McPhaden (2001), this bias in shortwave is probably due to using too high surface albedo over the ocean. In contrast, each component of the MERRA

flux, as well as the net flux, is within 10 W m⁻² of the TAO values in the western Pacific warm pool.

The tropical comparisons in Figs. 11 and 12 contrast with the global average results in Fig. 2 which show nearly all ORA-IP products gaining heat, probably at a higher rate than would be consistent with Argo data alone (e.g., Loeb et al. 2012).

## 5.2 Comparison with WHOI Stratus buoy

Outside the equatorial band a long-term (2000–2010) “Reference Timeseries” of in situ flux measurements from the WHOI Stratus buoy in the eastern subtropical South Pacific (mean location, 19.9°S, 85.3°W), has recently been made available to us (Bob Weller personal communication, 2014). The computation of the Stratus buoy fluxes uses hourly near-surface meteorological data and the COARE3.0b flux algorithm (Fairall et al. 1996). The height of the sensors above the sea surface is 2.3 m for humidity
<figure><img src="image_30.png" /><figcaption>**Fig. 13** Monthly timeseries of observed surface heat fluxes at the Stratus buoy (dashed blue lines) and the corresponding estimates derived from the ensemble of 11 ORA-IP products with all component fluxes available, see Table 1 (red lines with the pink shading</figcaption></figure>

indicating ensemble spread). Also shown to the right are the mean seasonal cycles in heat flux components over the period 2001–2009. Positive represents heat flux into the ocean. Units are in W m⁻²

and air-temperature, which compares favourably with the 2 m height at which temperature and humidity data are available from reanalyses. In addition, stability-dependent height corrections are applied to the buoy winds in order to adjust the 2.7 m anemometer measurements to 10 m height winds from reanalyses. Details of the buoy instrumentation and the bulk flux algorithm, along with the accuracy of the

computed buoy fluxes, are documented in Weller (2015) (see also http://uop.whoi.edu/projects/Stratus/stratus.html). In this comparison we use monthly mean fluxes from the Stratus buoy over a 9-year measurement period from January 2001 through December 2009 inclusive (108 months), within the ORA-IP time frame 1993–2009.
**Fig. 14** Monthly differences from the Stratus buoy, Ensemble ORA-IP (ERAi)-Stratus Buoy, using only the six ORA-IP products forced with ERAi fields using the CORE bulk formula: a radiative flux components and net radiation (short plus long wave heat fluxes), and b turbulent (sensible plus latent) and net heat fluxes along with surface wind speeds. Positive differences in surface wind speed (in m s⁻¹), ERAi-Buoy, indicate overestimation of ERAi compared to the buoy data. Positive flux differences (in W m⁻²) indicate larger ocean heat gain (shortwave and net fluxes) in the ORA-IP products and negative differences indicate larger ocean heat losses (longwave, latent and sensible fluxes) in ORA-IP compared to the buoy values

<figure><img src="image_31.png" /></figure>

Figure 13 shows monthly timeseries of observed surface heat fluxes at the Stratus buoy (dashed blue lines) and the corresponding estimates derived from the ensemble of 11 ORA-IP products extracted from the grid points nearest to the buoy. The Stratus 9-year mean net heat flux into the ocean is 38 W m⁻², with shortwave heating (191 W m⁻²) balanced by cooling dominated by latent heat (−103 W m⁻²), and to a lesser extent longwave radiation (−43 W m⁻²). The accuracy of the monthly mean surface heat fluxes from the Stratus buoy data are estimated by Weller (2015) to be: ±3 W m⁻² for net shortwave heat flux; ±2 W m⁻² for net longwave heat flux; ±4 W m⁻² for latent heat flux; and ±1.5 W m⁻² for sensible heat flux. There are compensating errors in flux components and the error in net surface heat flux is estimated to be ±8 W m⁻².

The ORA-IP ensemble reproduces the seasonal cycle in net heat flux dominated by variability in net shortwave radiation (Fig. 13a, b), although the net fluxes in the models are systematically lower than the buoy, with ~16 W m⁻² (42 %) less heat going into the ocean, nearly all of which (14 W m⁻²) is additional latent heat cooling. This clearly exceeds the error bounds of the computed net heat fluxes from the Stratus buoy by at least 50 %. Figure 13b–e show that seasonal component differences among ORA-IP products can be large: shortwave flux is overestimated, especially in spring and summer (October–January; Fig. 13b) and is compensated by more cooling due to latent, longwave and sensible fluxes in all seasons (Fig. 13c–e).

Figure 14a, b show ORA-IP ensemble differences from the buoy flux components using only the six ORA-IP products forced by ERAi using CORE bulk formula. Monthly differences in shortwave radiation are out of phase with longwave radiation (anti-correlation of −0.67), so that the net radiation bias is reduced in Fig. 14a. Differences in net flux are strongly correlated (up to 0.9) with differences in latent heat flux which dominate the turbulent fluxes, indicating that latent heat fluxes explain nearly all the variability in the differences between the buoy and the ORA-IP ensemble. The latent heat flux differences are in turn primarily explained by differences in surface wind speed (grey curve in Fig. 14b), with the ERAi winds being systematically too high (positive differences) compared to the buoy data, and with flux difference having significant correlated high frequency variability (~−0.6) associated with wind speed differences between ERAi and the buoy measurements.

All other bulk variables, including SST, air temperature and relative humidity (not shown) show positive biases (i.e., overestimation compared to the Stratus buoy observations) during 2001–2009. Mean offsets for SSTs are between 0.7 and 1.3 °C, but these differences only weakly correlate (<0.2) with latent flux biases in the ORA-IP ensemble, indicating that SSTs make only a small contribution to the differences in latent and net heat fluxes in the ocean reanalysis estimates.

Coincident with the variability discussed above, Weller (2015) found that increases in Stratus buoy wind speeds
and latent heat fluxes over the period 2001–2009 (primarily in spring and fall) were accompanied by a decrease in net heat flux of 39 W m⁻² or 104 % of the 9-year mean. A more extensive comparison with the WHOI Stratus buoy is now underway, trying to identify interannual signals more clearly, including looking at linear trends in ORA-IP and other gridded surface flux products over 2001–2009 for comparison with the buoy observations.

# 6 Summary and discussion

This paper looks at the surface heat fluxes from ocean reanalyses taking part in the GODAE/CLIVAR-GSOP Intercomparison Project ORA-IP. The emphasis is on the degree of agreement between products and, therefore, ensemble results are shown for means and spreads over the period 1993–2009, along with signal to noise ratio results to assess flux variability. Comparisons with available surface heat flux datasets, including local buoy measurements, are also made to highlight differences in input data and methodology.

Global mean heat fluxes (Fig. 2) generally show small net positive bias (i.e., heat flux into the ocean) compared with larger biases in observational based products (e.g., NOC2.0 and ISCCP/OAFlux). For products where data are available, the positive flux bias is largely compensated by the assimilation increments removing heat globally. Total heat gains (surface flux + assimilation increments) are typically 1–2 W m⁻², much smaller than most atmospheric reanalyses, but still larger than is realistic from Earth’s energy budget considerations (Loeb et al. 2012). We also show that compensation between resolved turbulent (latent and sensible) heat fluxes and assimilation heat loss largely takes place in the upper 100 m (Fig. 3).

The global compensation between surface heat fluxes and assimilation increments in near surface layers suggests that assimilation increments do contain information about surface heat flux errors that might be used in a similar way to Stammer et al. (2004), who use 4DVAR assimilation to generate surface heat flux corrections. Indeed, the GECCO2 product shown here already contains such corrections leading to a nearly closed global budget in Fig. 2. However, using the assimilation increments regionally and in time for flux correction would present many challenges. The increments are also correcting for errors in both horizontal and vertical heat transports (mixing errors near the surface), and will also depend on the distribution of available data to assimilate. An example of such transport errors shows up at the Equator in several products. In general the implied meridional heat transports are brought into much closer agreement with independent observations by combining the surface fluxes and assimilation increments

(Fig. 5b); however, large anomalous transports around the equator in some products can be seen where assimilation is unable to correct thermocline depths, suggesting that better bias corrections, e.g., Bell et al. (2004), may be needed.

Variability amongst ensemble members exhibit anticorrelations between short and long wave flux components, and between shortwave and latent heat flux components, consistent with these terms dominating local heat budgets nearly everywhere (i.e., ocean heat flux divergence variations are smaller). The ensemble mean seasonal cycle is highly consistent between members, with most areas showing monthly noise spread <10 W m⁻² (Fig. 6). However, the interannual signal/noise ratios in net surface fluxes are generally less than one in most areas, except in the tropical Pacific, where ENSO introduces large interannual signals which are captured well in the ensemble, Fig. 7e, f. We also show consistency in surface flux anomalies between ORA-IP and satellite radiation combined with OAFlux turbulent fluxes during two other interannual events, a large negative PDO in the Pacific in 2008, and a persistent negative NAO period in the Atlantic in the later part of 2009 (Figs. 9, 10).

The reanalysis products are also compared to tropical buoy fluxes for 2007–2009 (when all flux components from eight tropical buoys are available), and to a long-term (2001–2009) reference dataset of fluxes at the Stratus buoy in the eastern south Pacific. All the reanalysis products show underestimation of ocean heat gain at the tropical buoys (one-third smaller Qnet than observed) primarily due to latent heat flux errors and, to a lesser degree, short and longwave radiation errors (Fig. 11). At the Stratus buoy, the ORA-IP ensemble again underestimates ocean heat gain and the temporal variability of this bias is completely dominated by latent heat flux errors caused by differences in the surface winds imposed in the reanalysis models compared with those measured directly at the buoy (Fig. 14b). This is a strong indication that better surface winds are a likely prerequisite for better surface fluxes. Given the strong relationships between SST gradients and surface winds this suggests coupled reanalyses may provide improved results (see further below).

These reanalysis products were not designed with air-sea fluxes in mind and the assessment here is relying on the influence of upper ocean temperatures, which are being controlled through data assimilation, to positively influence the surface heat fluxes. Song and Yu (2013) uses a cage or pool-area heat budget analysis together with in situ surface and subsurface measurements to examine the consistency of nine flux climatologies (from atmospheric reanalyses or blended products) and to identify uncertainties, in the western Pacific warm pool region. Such diagnostics using in situ ocean data have some potential to validate surface flux products and remove flux biases, but can only be used in regions with small lateral exchanges (e.g., the Pacific warm pool), if these lateral exchanges are assumed to be unknown.
Other efforts to improve air-sea fluxes are documented in Trenberth et al. (2011), Mayer et al. (2014) or Liu et al. (2015), who use atmospheric reanalysis products to assess air-sea fluxes based on atmospheric transport divergences. This approach could also work in the ocean, since the ocean currents should be strongly geostrophically constrained by the assimilated Argo profile data, although the ability of the ocean to store heat means that, to 1st order, surface flux errors are more directly compensated by assimilation increments, as shown here in Fig. 3. The comparison of ocean heat transports is a future objective of the ORA-IP program and we would argue that a demonstrable ability to reproduce ocean heat and freshwater transports with smaller error bounds is a key requirement of a good ocean assimilation product since this is information that is not directly available from the upper ocean observations of temperature and salinity.

In future coupled atmosphere-ocean reanalyses using coupled data assimilation schemes should help improve air-sea fluxes within reanalysis products. Truly coupled approaches to reanalysis would allow critical observational data such as sea surface temperatures to be assimilated correctly, e.g., from the ESA CCI program, Merchant et al. (2014), including their error representations. In the current ORA-IP products there are several coupled reanalysis results, but lower resolution models are being used and SSTs are generally still being assimilated as statistically gridded surface products, which will not correctly represent the spatial and temporal information content available from the satellite observational fields.

Finally, we note that the availability of in situ flux datasets suitable for evaluating large scale long-term products such as from reanalyses is still very limited (Yu et al. 2013). The tropical buoys used here provide a suitable first step but more mid-latitude buoys, e.g., from OceanSITES, with longer time records are needed in order to assess the lower frequency variability in heat fluxes outside the tropics. Cross calibration between buoy deployments is needed to create these long record datasets, as provided for the Stratus buoy (Bob Weller personal communication), and the wider availability of such products would allow greater use by the modelling and satellite data processing communities.

**Acknowledgments** We would like to thank Bob Weller for making available the WHOI Stratus buoy dataset, and Simon Josey and Pat Hyder for useful discussions on air-sea fluxes. Comments and suggestions provided by Simon Josey and two anonymous reviewers contributed to an improved presentation of these results. This work was supported by the Natural Environment Research Council (NERC) under the National Centre for Earth Observation (NCEO) program.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## References

Andersson A, Fennig K, Klepp C, Bakan S, Grassi H, Schulz J (2010) The Hamburg ocean atmosphere parameters and fluxes from satellite data—HOAPS-3. Earth Syst Sci Data 2:215–234. doi:10.5194/essd-2-215-2010

Arndt DS, Baringer MO, Johnson MR (2010) State of the climate in 2009. Bull Am Meteor Soc 91:s1–s222. doi:10.1175/BAMS-91-7-StateoftheClimate

Bacon S (1997) Circulation and fluxes in the North Atlantic between Greenland and Ireland. J Phys Oceanogr 27(7):1420–1435

Balmaseda MA, Dee D, Vidard A, Anderson DLT (2007) A multivariate treatment of bias for sequential data assimilation: application to the tropical oceans. Q J R Meteorol Soc 133:167–179

Balmaseda MA, Mogensen K, Weaver AT (2013a) Evaluation of the ECMWF ocean reanalysis system ORAS4. Q J R Meteorol Soc 139:1132–1161. doi:10.1002/qj.2063

Balmaseda MA, Trenberth KE, Källén E (2013b) Distinctive climate signals in reanalysis of global ocean heat content. Geophys Res Lett 40:1754–1759. doi:10.1002/grl.50382

Balmaseda M et al (2015) The Ocean Reanalyses Intercomparison Project (ORA-IP). J Oper Oceanogr 7:81–99

Behringer D (2007) The global ocean data assimilation system at NCEP. In: 11th symposium on integrated observing and assimilation systems for atmosphere, oceans and land surface. American Meteorological Society

Bell MJ, Martin MJ, Nichols NK (2004) Assimilation of data into an ocean model with systematic errors near the equator. Q J R Meteorol Soc 130:873–893. doi:10.1256/qj.02.109

Bentamy A, Grodsky SA, Katsaros K, Mestas-Nunez AM, Blanke B, Desbiolles F (2013) Improvement in airsea flux estimates derived from satellite observations. Int J Remote Sens 34(14):5243–5526

Berry DI, Kent EC (2009) A new air-sea interaction gridded dataset from ICOADS with uncertainty estimates. Bull Am Meteorol Soc 90(5):645–656. doi:10.1175/2008BAMS2639.1

Blockley EW, Martin MJ, McLaren AJ, Ryan AG, Waters J, Lea DJ, Mirouze I, Peterson KA, Sellar A, Storkey D (2014) Recent development of the Met Office operational ocean forecasting system: an overview and assessment of the new global FOAM forecasts. Geosci Model Dev Discuss 6:6219–6278. doi:10.5194/gmdd-6-6219-2013

Bourlès B et al (2008) The PIRATA program: history, accomplishments, and future directions. Bull Am Meteorol Soc 89:1111–1125

Bryden HL, Imawaki S (2001) Ocean heat transport. In: Siedler G, Church J, Gould J (eds) Ocean circulation and climate: observing and modelling the global ocean. Academic Press, San Francisco, pp 455–474 (International Geophysics Series 77)

Cayan DR (1992) Latent and sensible heat flux anomalies over the northern oceans: driving the sea surface temperature. J Phys Oceanogr 22(8):859–881

Chang Y-S, Zhang S, Rosati A, Delworth TL, Stern WF (2013) An assessment of oceanic variability for 1960–2010 from the GFDL ensemble coupled data assimilation. Clim Dyn 40:775–803. doi:10.1007/s00382-012-1412-2

Dee DP et al (2011) The ERA-Interim reanalysis: configuration and performance of the data assimilation system. Q J R Meteorol Soc 137:553–597. doi:10.1002/qj.828

Dussin R, Barnier B (2013) The making of DFS 5.1. DRAKKAR Report, February 2013

Fairall CW, Bradley EF, Rogers DP, Edson JB, Young GS (1996) Bulk parametrization of air-sea fluxes for tropical ocean global atmosphere coupled ocean-atmosphere response experiment. J Geophys Res 101:3747–3764. doi:10.1029/95JC03205
Fairall CW, Bradley EF, Hare JE, Grachev AA, Edson JB (2003) Bulk parameterization of air-sea fluxes: updates and verification for the COARE algorithm. J Clim 16(4):571–591

Ferry N, Parent L, Garric G, Bricaud C, Testut C-E, Le Galloude O, Lellouche J-M, Drevillon M, Greiner E, Barnier B, Molines J-M, Jourdain NC, Guinehut S, Cabanes C, Zawadzki L (2012) GLO-RYS2V1 global ocean reanalysis of the altimetric era (1992–2009) at meso scale. Mercator Q Newsl 44:29–39. http://www.mercator-ocean.fr/eng/actualites-agenda/newsletter/newsletter-Newsletter-44-Various-areas-of-benefit-using-the-Mercator-Ocean-products. GLORYS2V1_MercatorNewsletter44.pdf

Fujii Y, Nakaegawa T, Matsumoto S, Yasuda T, Yamanaka G, Kamachi M (2009) Coupled climate simulation by constraining ocean fields in a coupled model with ocean data. J Climate 22:5541–5557. doi:10.1175/2009JCLI2814.1

Fujii Y, Tsujino H, Toyoda T, Nakano H (2015) Enhancement of the southward return flow of the Atlantic Meridional Overturning Circulation by data assimilation and its influence in an assimilative ocean simulation forced by CORE-II atmospheric forcing. Clim Dyn (Submitted to the same special issue)

Ganachaud A, Wunsch C (2003) Large-scale ocean heat and freshwater transports during the World Ocean Circulation Experiment. J Clim 16(4):696–705

Griffies SM et al (2009) Coordinated ocean-ice reference experiments (COREs). Ocean Model 26(1):1–46

Haines K, Valdivieso M, Zuo H, Stepanov VN (2012) Transports and budgets in a 1/4 global ocean reanalysis 1989–2010. Ocean Sci Discuss 9(1):261–290

Josey SA, Smith SR (2006) Guidelines for evaluation of air-sea heat, freshwater and momentum flux datasets. CLIVAR Global Synthesis and Observations Panel (GSOP) Report 12

Josey SA, Kent EC, Taylor PK (1999) New insights into the ocean heat budget closure problem from analysis of the SOC air-sea flux climatology. J Clim 12(9):2856–2880

Josey SA, Gulev S, Yu L (2013) Exchanges through the ocean surface. In: Siedler G, Griffies S, Gould J, Church J (eds) Ocean circulation and climate, 2nd edn. A 21st century perspective. Academic Press, International Geophysics Series, vol 103, pp 115–140

Josey SA, Lisan Y, Gulev S, Xiangze Jin J, Tilinina N, Barnier B, Brodeau L (2014) Unexpected impacts of the Tropical Pacific array on reanalysis surface meteorology and heat fluxes. Geophys Res Lett 41(17):6213–6220

Kallberg P (1987) Forecast drift in ERA-Interim. Technical Report ERA Report Series, No 10, ECMWF

Kalanay E et al (1996) The NCEP/NCAR 40-year reanalysis project. Bull Am Meteorol Soc 77:437–472

Kanamitsu M, Ebitsuzaki W, Woolen J, Yang SK, Hnilo JJ, Fiorino M, Potter G (2002) NCEP-DOE AMIP-II reanalysis (R-2). Bull Am Meteorol Soc 83:1631–1643

Kobayashi S et al (2015) The JRA-55 reanalysis: general specifications and basic characteristics. Accepted for publication in Journal of the Meteorological Society of Japan (JMSJ) 93(1). doi:10.2151/jmsj.2015.001

Köhl A (2015) Evaluation of the GECCO2 ocean synthesis: transports of volume, heat and freshwater in the Atlantic. Q J R Meteorol Soc 141:166–181. doi:10.1002/qj.2347

Kubota M, Iwasaka N, Kizu S, Konda M, Kutsuwada K (2002) Japanese ocean flux data sets with use of remote sensing observations (J-OFURO). J Oceanogr 58(1):213–225

Large WG, Pond S (1981) Open ocean momentum flux measurements in moderate to strong winds. J Phys Oceanogr 11:324–336

Large WG, Yeager SG (2004) Diurnal to decadal global forcing for ocean and sea-ice models: the data sets and flux climatologies. Technical Report, TN-460 + STR, NCAR, 105 pp

Large WG, Yeager SG (2009) The global climatology of an interannually varying air-sea flux data set. Clim Dyn 33(2–3):341–364

Lellouche J-M et al (2013) Evaluation of global monitoring and forecasting systems at Mercator Ocean. Ocean Sci 9:57–81. doi:10.5194/os-9-57-2013

Liu C, Allan RP, Berrisford P, Mayer M, Hyder P, Loeb N, Smith D, Vidale P-L, Edwards JM (2015) Combining satellite observations and reanalysis energy transports to estimate global net surface energy fluxes 1985–2012. Accepted in J Geophys Res. doi:10.1002/2015JD023264/pdf

Loeb NG, Wielicki BA, Doelling DR, Smith GL, Keyes DF, Kato S, Manalo-Smith N, Wong T (2009) Toward optimal closure of the Earth’s top-of-atmosphere radiation budget. J Clim 22(3):748–766. doi:10.1175/2008JCLI2637.1

Loeb NG, Lyman JM, Johnson GC, Allan RP, Doelling DR, Wong T, Soden BJ, Stephens GL (2012) Observed changes in top-of-the-atmosphere radiation and upper-ocean heating consistent within uncertainty. Nat Geosci 5(2):110–113. doi:10.1038/ngeo1375

Lumpkin R, Speer K (2007) Global ocean meridional overturning. J Phys Oceanogr 37(10):2550–2562

Macdonald A, Baringer M (2013) Ocean heat transport. In: Siedler G, Griffies SM, Gould J, Church JA (eds) Ocean circulation and climate: a 21st century perspective. International geophysics series, vol 103. Academic Press, pp 759–785

Mayer M, Haimberger L, Balmaseda MA (2014) On the energy exchange between tropical ocean basins related to ENSO. J Clim 27:6393–6403. doi:10.1175/JCLI-D-14-00123.1

McPhaden M, Busalacchi A, Cheney R, Donguy J-R, Gage K, Halpern D, Ji M, Meyers PJG, Mitchum G, Niiler P, Picaut J, Reynolds R, Smith N, Takeuchi K (1998) The tropical ocean global atmosphere observing system: a decade of progress. J Geophys Res 103:14169–14240

McPhaden M et al (2009) RAMA: the research moored array for African-Asian-Australian monsoon analysis and prediction. Bull Am Meteorol Soc 90:459–480

Merchant CJ, Embury O, Roberts-Jones J, Fiedler E, Bulgin CE, Corlett GK, Good S, McLaren A, Rayner N, Morak-Bozzo S, Donlon C (2014) Sea surface temperature datasets for climate applications from Phase 1 of the European Space Agency Climate Change Initiative (SST CCI). Geosci Data J. doi:10.1002/gdj.3.20

Palmer MD, McNeall DJ (2014) Internal variability of Earth’s energy budget simulated by CMIP5 climate models. Environ Res Lett 9:034016. doi:10.1088/1748-9326/9/3/034016

Peixoto JP, Oort AH (1992) Physics of climate. American Institute of Physics, New York, p 520

Peterson TC, Baringer MO (2009) State of the climate in 2008. Bull Am Meteorol Soc 90:S1–S196. doi:10.1175/BAMS-90-8-StateoftheClimate

Rienecker MM, Suarez MJ, Gelaro R, Todling R, Bacmeister J, Liu E, Bosilovich MG et al (2011) MERRA: NASA’s modern-era retrospective analysis for research and applications. J Clim 24(14):3624–3648

Roemmich D et al (2015) Unabated planetary warming and its ocean structure since 2006. Nat Clim Change. doi:10.1038/nclimate2514

Schiller A, Lee T, Masuda S (2013) Methods and applications of ocean synthesis in climate research. In: Siedler G, Griffies SM, Gould J, Church JA (eds) Ocean circulation and climate: a 21st century perspective. International geophysics series, vol 103. Academic Press, pp 581–601

Servain J, Busalacchi AJ, McPhaden MJ, Moura AD, Reverdin G, Vianna M, Zebiak SE (1998) A pilot research moored array in the tropical Atlantic (PIRATA). Bull Am Meteorol Soc 79(10):2019–2031

Song X, Yu L (2013) How much net surface heat flux should go into the Western Pacific Warm Pool? J Geophys Res Oceans 118:3569–3585. doi:10.1002/jgrc.20246
Stammer D, Ueyoshi K, Köhl A, Large WG, Josey SA, Wunsch C (2004) Estimating air–sea fluxes of heat, freshwater, and momentum through global ocean data assimilation. J Geophys Res 109:C05023. doi:10.1029/2003JC002082

Storto A, Masina S, Dobricic S (2014) Estimation and impact of nonuniform horizontal correlation length scales for global ocean physical analyses. J Atmos Ocean Technol 31:2330–2349

Storto et al (2015) Steric sea level variability (1993–2010) in an ensemble of ocean reanalyses and objective analyses. Clim Dyn. Early Online Release. doi:10.1007/s00382-015-2554-9

Toyoda T, Fujii Y, Yasuda T, Usui N, Iwao T, Kuragano T, Kamachi M (2013) Improved analysis of the seasonal-interannual fields by a global ocean data assimilation system. Theor Appl Mech Jpn 61:31–48

Trenberth KE, Fasullo JT, Mackaro J (2011) Atmospheric moisture transports from ocean to land and global energy flows in reanalyses. J Clim 24:4907–4924. doi:10.1175/2011JCLI4171.1

Valdivieso M, Haines K, Zuo H, Lea D (2014a) Freshwater and heat transports from global ocean synthesis. J Geophys Res Oceans 119(1):394–409. doi:10.1002/2013JC009357

Valdivieso M, Haines K, Balmaseda M, Barnier B, Chang Y, Ferry N, Fujii Y, Köhl A, Lee T, Martin M, Storto A, Toyoda T, Wang X, Waters J, Xue Y, Yin Y (2014) Heat fluxes from ocean and coupled reanalyses. CLIVAR exchanges, special issue: ongoing efforts on ocean reanalyses intercomparison, no. 64, vol 19(1), pp 28–31

von Schuckmann et al (2015) Consistency between planetary energy balance and ocean heat storage (CONCEPT-HEAT). A prospec- tus for the CLIVAR research focus. WCRP Report No. 14/2015, CLIVAR Report No. 203, June 2015

Wang W, McPhaden MJ (2001) What is the mean seasonal cycle of surface heat flux in the equatorial Pacific? J Geophys Res Oceans (1978–2012) 106(C1):837–857

WCRP (2012) Action plan for WCRP research activities on surface fluxes. WCRP Informal/Series Report No. 01, January 2012

Weller RA (2015) Variability and trends in surface meteorology and air–sea fluxes at a site off Northern Chile. J Clim. doi:10.1175/ JCLI-D-14-00591.1

WGASF (2000) Intercomparison and validation of ocean–atmosphere energy flux fields. In: Taylor PK (ed) Final report of the Joint WCRP/SCOR Working Group on Air–Sea Fluxes (WGASF), WCRP-112, WMO/TD-1036, 306 pp

Wunsch C, Heimbach P (2013) Ocean circulation and climate—a 21st century perspective. In: Sielder G, Church J, Gri.Jes S, Gould J, Church J (eds) Dynamically and kinematically consistent global ocean circulation and ice state estimates. International geophys- ics series, vol 103. Academic Press, Elsevier

Xue Y, Huang B, Hu Z-Z, Kumar A, Wen C, Behringer D, Nadiga S (2011) An assessment of oceanic variability in the NCEP climate forecast system reanalysis. Clim Dyn 37(11–12):2511–2539

Xue Y et al (2012) A comparative analysis of upper-ocean heat con- tent variability from an ensemble of operational ocean reanaly- ses. J Clim 25(20):6905–6929

Yin Y, Alves O, Oke P (2011) An ensemble ocean data assimilation system for seasonal prediction. Mon Wea Rev 139:786–808

Yu L, Jin X, Weller RA (2008) Multidecade global flux datasets from the objectively analyzed air–sea fluxes (OAFlux) project: latent and sensible heat fluxes, ocean evaporation, and related surface meteorological variables. Technical Report OAFlux Project (OA2008–01), Woods Hole Oceanographic Institution

Yu L, Haines K, Bourassa M, Cronin M, Gulev S, Josey S, Kato S, Kumar A, Lee T, Roemmich D (2013) Towards achieving global closure of ocean heat and freshwater budgets: Recommendations for advancing research in air–sea fluxes through collaborative activities. International CLIVAR Project Office, 2013: Interna- tional CLIVAR Publication Series No 189. http://www.clivar. org/sites/default/files/ICPO189_WHOI_fluxes_workshop.pdf

Zhang Y, Rossow WB, Lacis AA, Oinas V, Mishchenko MI (2004) Calculation of radiative fluxes from the surface to top of atmos- phere based on ISCCP and other global data sets: refinements of the radiative transfer model and the input data. J Geophys Res 109:D19105. doi:10.1029/2003JD004457