Inverse Modeling of Global and Regional Energy and Water Cycle Fluxes using Earth Observation Data

The NASA Energy and Water Cycle Study (NEWS) climatology is a self-consistent coupled annual and seasonal cycle solution for radiative, turbulent, and water ﬂuxes over Earth’s surface using Earth observation data covering 2000–09. Here we seek to improve the NEWS solution, particularly over the ocean basins, by considering spatial covariances in the observation errors (some evidence for which is found by comparing ﬁve turbulent ﬂux products over the oceans) and by introducing additional horizontal transports from ocean reanalyses as weak constraints. By explicitly representing large error covariances between surface heat ﬂux components over the major ocean basins we retain the ﬂux contrasts present in the original data and infer additionalheatlossesovertheNorthAtlanticOcean,moreconsistentwitha strongAtlanticoverturning.This change does not alter the global ﬂux balance but if only the errors in evaporation and precipitation are correlated then those ﬂuxes experience larger adjustments (e.g., the surface latent heat ﬂux increases to 85 6 2 Wm 2 2 ). Replacing SeaFlux v1 with J-OFURO v3 (Japanese Ocean Flux Data Sets with Use of Remote Sensing Observations) ocean ﬂuxes also leads to a considerable increase in the global latent heat loss as well as a larger North Atlantic heat loss. Furthermore, including a weak constraint on the horizontal trans-ports of heat and freshwater from high-resolution ocean reanalyses improves the net ﬂuxes over the North Atlantic,CaribbeanSea,andArcticOcean,withoutanyimpactontheglobalﬂuxbalances.Theseresultssuggestthat better characterized ﬂux uncertainties can greatly improve the quality of the optimized ﬂux solution.


Introduction
Vertical and horizontal energy flows between Earth's surface, atmosphere, and space play a fundamental role in establishing the large-scale atmosphere and ocean circulation patterns driving weather and climate. The water cycle is closely coupled to these energy flows due to the exchanges of latent heat that occur during evaporation/transpiration and precipitation.
A wide variety of Earth observation (EO) datasets are now available, enabling different vertical components of the Earth's energy and water cycles to be quantified. For example, spatial and temporal variations in the top-ofatmosphere (TOA) radiation fluxes, both shortwave and longwave, are now well monitored by satellites, including the CERES product (Wielicki et al. 1996;Stephens et al. 2012). The vertical structure of clouds derived from CloudSat/CALIPSO observations (Henderson et al. 2013) permits surface radiative fluxes to be estimated reasonably precisely. Ocean heat content estimates from Argo measurements (Roemmich et al. 2015) allow longterm oceanic energy storage to be included in a determination of the global energy cycle. Satellite-based turbulent flux products can be used to determine surface latent and sensible heat exchanges. Precipitation fluxes are available from combinations of satellite observations such as that produced in the GPCP project (Adler et al. 2003;Huffman et al. 2009). Finally there are over a decade of surface water storage estimates from the GRACE satellites (Tapley et al. 2004;Chambers and Bonin 2012;Johnson and Chambers 2013).
Large imbalances have been shown to exist when independent energy flux observations are combined at Denotes content that is immediately available upon publication as open access. both global and regional scales (e.g., Josey et al. 1999); the net vertical fluxes are inconsistent with any realistic storage or horizontal transports, limiting the value of such observations for constraining climate models (Wild et al. 2015). Attempts to reconcile these imbalances have had rather limited impact (e.g., Grist and Josey 2003). An alternative approach that uses reanalysis output to estimate energy and water fluxes, particularly at the surface, has also been widely employed (e.g., Trenberth and Caron 2001;Fasullo and Trenberth 2008a,b;Liu et al. 2015Liu et al. , 2017Valdivieso et al. 2017). However, estimates of fluxes based on these reanalyses alone will be sensitive to numerical model parameterizations and biases, and subtle water mass corrections should be applied prior to the energy flux calculations (Mayer et al. 2017). Furthermore, these estimates will generally lack robust uncertainty estimates as most current operational centers do not provide such uncertainties, although ensemble and multiproduct approaches partly mitigate this concern.
A novel approach to reconciling the satellite derived observational fluxes was developed for the NASA Energy and Water Cycle Study (NEWS). Whereas EO datasets have previously been used for assessing either energy (Isemer et al. 1989;da Silva et al. 1994;Grist and Josey 2003) or water (Trenberth et al. 2011) budgets separately, the NEWS papers by L' Ecuyer et al. (2015) and Rodell et al. (2015b) demonstrate for the first time how to explicitly couple the energy and water cycles. The NEWS method takes the 2000-09 EO-based estimates of vertical fluxes and their uncertainties, at both the top of the atmosphere (TOA) and Earth's surface, in 16 land and ocean regions across the globe. A variety of physical balance constraints are then imposed and the fluxes are allowed to change according to their uncertainties in order to close the energy and water budgets in the atmosphere and at the surface. Both annual and monthly 10-yr mean solutions are determined.
The main physical balance constraints used in the NEWS solution are as follows. Land areas are assumed not to absorb any energy in the long-term mean, and the global ocean is assumed to absorb ;1.0 6 0.6 W m 22 in line with Argo measurements, although the spatial distribution of this warming is not imposed. The atmosphere is assumed to absorb very little energy due to its low heat capacity. Water is conserved within the Earth system, providing additional constraints. An estimate of surface (and groundwater) runoff from the land regions into ocean basins is incorporated into the fit, and this is complemented by horizontal water exchange between regions in the atmosphere provided by the MERRA reanalysis (Rienecker et al. 2011). These are the only horizontal exchanges that have specified prior values; no estimates of horizontal energy exchange are employed. The optimized solution then provides estimates of all of the energy and water exchanges within the atmosphere and oceans, and at Earth's surface, along with new uncertainty estimates.
The purpose of the current paper is to extend the NEWS approach in two ways. The first is to explicitly permit spatially correlated errors in the original EO flux products. We show that this can make the NEWS solutions more realistic, especially over the oceans, bringing the regional surface fluxes closer to reanalysis-derived values and modifying the global estimates of turbulent heat exchanges. Second we bring reanalysis results into the NEWS framework by using ocean reanalysis-derived horizontal energy and water flux priors to complement the vertical fluxes taken from EO data. Such a complementary combination of approaches is advocated, for example, by the recent CLIVAR Research Focus on Planetary heat balance and ocean heat storage (CONCEPT-HEAT) international program looking at Earth's energy imbalance.
The paper is organized as follows. Section 2 describes the original NEWS flux solution, including the key observation datasets and the solution method, and compares surface fluxes and transports between the NEWS solution and reanalysis-derived quantities. Section 3 details the introduction of explicit error covariances over the oceans and the corresponding impact on the regional and global fluxes. We demonstrate the evidence for error covariances in five different ocean turbulent heat flux products and also explore the impact of using different turbulent fluxes on the optimized solution. In section 4 the incorporation of ocean reanalysis values into the solution is described and its influence on the regional and global fluxes is evaluated. Section 5 then provides some discussion of the new optimized flux fields and their uncertainties, as well as of additional improvements that could be made. Finally section 6 contains a summary and conclusions.

a. Energy and water flux input data
We use very similar flux datasets to those employed in the NEWS study, which were chosen based on a variety of factors such as the time period covered, the availability of robust uncertainty estimates, and the extent to which the data had been adopted by the scientific community. The datasets are extensively discussed in L' Ecuyer et al. (2015) and Rodell et al. (2015b) and published online at Rodell et al. (2015a), and so will be only briefly reviewed here. The data cover the period 2000-09 (inclusive) and both the 10-yr annual mean and 10-yr mean seasonal cycle data are used. We also use the uncertainties provided for both time scales. The time period chosen has several desirable properties: for example, a large variety of datasets, both satellite and otherwise, are available for uncertainty evaluation, and the trend in global mean surface temperature was small (although other aspects of climate change were evident during this time). Although some datasets have missing years, the interannual variations are smaller than the uncertainties so the decadal mean can be taken as representative of the entire period.
TOA and surface radiative fluxes are taken from the CERES dataset (Wielicki et al. 1996) and precipitation measurements are obtained from the Global Precipitation Climatology Product v2.2 (Adler et al. 2003;Huffman et al. 2009). Ocean evaporation is taken from the SeaFlux v1.0 product (Curry et al. 2004), and land evapotranspiration is taken from a dataset produced at Princeton University (Vinukollu et al. 2011) combined with the MERRA (Rienecker et al. 2011;Bosilovich et al. 2011) and GLDAS (Rodell et al. 2004) reanalyses produced by NASA. Observations of runoff from land are obtained from a dataset produced at the University of Washington (Clark et al. 2015). Terrestrial water storage changes are calculated from globally detrended GRACE satellite measurements (Tapley et al. 2004;Chambers and Bonin 2012;Johnson and Chambers 2013). Finally, total precipitable water vapor is calculated from the AIRS and AMSR-E instruments on the Aqua satellite (Fetzer et al. 2006).
We favor solely using the direct CERES observations of TOA radiation but the NEWS team also employed the Global Energy and Water Cycle Experiment Surface Radiation Budget dataset (Gupta et al. 1999) and the International Satellite Cloud Climatology Project Flux Data product (Zhang et al. 2004). Furthermore we do not use priors for atmospheric water vapor convergence [NEWS took these from MERRA, QuikSCAT (Liu et al. 2006), and the Passive Microwave Water Cycle dataset (Hilburn 2009)]. This means we do not use water convergence data over individual ocean basins (Rodell et al. 2015b, their Table 3) until section 4 when we explicitly introduce additional ocean transport data. The motivation for initially excluding the reanalysis data is that the surface runoff observations already provide a strong enough constraint on the horizontal transport of water.
The 16 land and ocean regions used in our study are the same as those used in the NEWS solution, because we have used much of the same data provided by NEWS, and are presented in Fig. 1.

b. Inverse model
An inverse model is used to impose budget constraints on all the component fluxes of the annual mean seasonal cycles of energy and water in the 16 regions used in the NEWS study. The vertical and horizontal energy and water fluxes combined in the fit are as follows: ISR (incoming TOA radiation); OSR and OLR (outgoing shortwave/longwave TOA radiation); DSR and DLR (downwelling shortwave/longwave surface radiation); USW and ULW (upwelling shortwave/longwave surface radiation); SH and LE (upwelling surface sensible/latent heat flux); P (precipitation); ADC and AWC (atmospheric horizontal dry static energy/water convergence); SEC and SWC [(sub)surface horizontal energy/water convergence]; and dS and dW (change in subsurface/atmospheric water storage). The terms P and LE appear in both energy and water balance equations and ensure the energy and water cycles are coupled.
In each region r the following balance equations must be satisfied. The water cycle terms are converted to units of energy prior to including them in the equations.
The first balance equation expresses the fact that the annual average atmospheric energy residual is close to zero: ISR r 2 OLR r 2 OSR r 2 DLR r 2 DSR r 1 ULW r 1 USW r 1 SH r 1 P r 1 ADC 5 0 6 0:05 W m 22 "r . (1) The second equation represents (sub)surface energy balance: where the annual average b r 6 s r is 0 6 0.05 W m 22 and 1.0 6 0.6 W m 22 for land and oceans, respectively. The latter values reflect the global ocean energy imbalance estimated from Argo measurements (Roemmich et al. 2015). The two energy balance equations neglect the heat capacity of water being exchanged with the atmosphere; in other words, they do not account for evaporation and precipitation and horizontal water convergences taking place at different temperatures, but these terms are likely small compared to adjustments made to other fluxes (Mayer et al. 2017). On a monthly time scale there can be residual surface energy storage that is not distinguishable from changes in horizontal convergence (SEC r ) given the observations above. The equation expressing atmospheric water balance is LE r 2 P r 1 AWC 5 dW r "r and similarly the (sub)surface water balance is P r 2 LE r 1 SWC r 5 dS r "r .
On an annual average the storage change terms dS r and dW r are all assumed to be zero; in this case the variable AWC is redundant (becauseAWC r 5 SWC r " r ) and is thus not included in the fit. In our implementation of the fit the terms ADC, SEC, and AWC do not have prior estimates but are instead derived indirectly. The values of SWC and dS only have priors over land areas; section 4 describes the use of prior values of SWC (and SEC) over the oceans as additional constraints. Last, the global sum of each horizontal convergence is required to be exactly zero: å r ADC r 5 å r SEC r 5 å r AWC r 5 å r SWC r 5 0 (5) except when a convergence cannot be distinguished from a storage change, as occurs for SEC r in the monthly fit. A cost function is set up to determine the optimized fluxes, convergences, and storage terms. The vector F contains all optimized variables that have prior observed values (F obs ); the optimized variables are penalized in the cost function according to the offsets relative to their prior values weighted by the inverse of the observation error covariance matrix R obs . Variables without prior values are contained in the vector G and are determined by the various physical balance constraints. The cost function is where there are K physical budget constraints. Lagrange multipliers l are used to enforce any exact physical constraints. Inexact constraints are expressed by cost function terms such as where S(F, G) is a particular linear combination of optimized variables and V 6 s V is the value and uncertainty of the constraint on this combination [see Eqs.
(1) and (2)]. The minimum value of the cost function, J min , is found in order to determine the values of F and G that best match both F obs and the physical balances. A full description the cost function and its minimization can be found in appendix A. Uncertainties in the solution are determined by inverting the cost function Hessian matrix. These uncertainties contain contributions from the disparate datasets' uncertainties as well as the constraints that link the physical variables, and should not be interpreted as independent improvements in each observational flux uncertainty.
c. Net surface fluxes, ocean transports, and globally averaged fluxes There are some small differences in our implementation of this solution compared to the NEWS fit; as well as using slightly different datasets (section 2a), we perform the minimization directly in the 16 regions (assumed to be uncorrelated), in contrast to the NEWS approach of using a two-stage fit that is initially solved using a global ocean region. The regional fluxes determined in the inverse method are very similar to those shown in L' Ecuyer et al. (2015) and Rodell et al. (2015b) but with slightly smaller final uncertainties.
In this section we describe some of the problems with the NEWS solution, providing motivation for our work in the remainder of the paper. We first consider the time mean (2000-09) net downward surface energy flux (NSF, consisting of net downward longwave and shortwave radiation and sensible and latent heat fluxes). Figure 2 shows the NSF in the 16 NEWS regions both for the initial EO data (Fig. 2a) and from the optimized solution ( Fig. 2b).
The initial EO data clearly have large NSF imbalances over the land surface as well as a grossly positive energy flux into the oceans of approximately 20-30 W m 22 , which has been noted in many previous heat flux studies (e.g., Grist and Josey 2003). The inverse solution has successfully reduced the values of NSF to zero over the land areas, as required by Eq. (2), and reduced NSF over the ocean basins to much smaller values that are both positive and negative, consistent with a small global average. The uncertainties over land become small because of the strong NSF regional constraint but the uncertainties over the oceans have not greatly reduced, although the component fluxes (radiative and turbulent) that they consist of now have correlated uncertainties imposed by the optimization process.
For comparison, Fig. 3 shows the regional NSF based on the CERES TOA radiation fluxes combined with the ERA-Interim reanalysis (Dee et al. 2011) horizontal energy convergences, following the recent surface flux derivation approach of Liu et al. (2015). Only the NSF can be provided by this method, without any breakdown into flux components, and there is also no objective means of calculating an uncertainty in the final product. However, some idea of the uncertainties may be taken from the residual NSF into the land areas (about 64 W m 22 ), which presumably reflects inaccuracies in the imposed reanalysis transports. Liu et al. (2015) provide ad hoc solutions to these problems.
Leaving aside the land areas noted above, the reanalysis-derived NSF is clearly within the uncertainty limits of the NEWS solution over nearly all ocean basins, but large discrepancies can be seen over the North Atlantic, Caribbean, and Arctic basins, which exhibit large surface energy losses. The North Atlantic basin is known to be a region of strong heat loss from the oceans to the atmosphere due to the northward transport of heat in the Atlantic meridional overturning circulation (AMOC) (Hall and Bryden 1982;Ganachaud andWunsch 2000, 2001;Trenberth and Caron 2001;Mignac et al. 2018). In the NEWS solution the implied cross-equatorial heat flux is 0.01 6 0.53 PW southward and the heat flux into the Arctic is 0 6 130 TW (assuming no heat flux through the Bering Strait; see appendix B for a description of the method used to determine these quantities). The North Atlantic NSF in the NEWS solution and the associated cross-equatorial heat flux are clearly unrealistic. In contrast, the implied ocean cross-equatorial steady state heat flux in the Atlantic Ocean in the reanalysis-derived solution is 0.72 PW northward, consistent with a strong AMOC, and the heat flux into the Arctic is 120 TW. The crossequatorial transport in the reanalysis-derived solution is considerably more realistic, although the implied ocean heat transport into the Arctic is much too strong in the ERA-Interim solution and is rather more realistic in the NEWS solution. We will focus on the North Atlantic, Caribbean, and Arctic basins as we explore modifications to the NEWS solution process in order to provide improved objective flux analyses.
Considering the global fluxes in the NEWS solution, L'Ecuyer et al. (2015) comment that there are considerable variations between published estimates of some of the components. The global sensible heat flux mean of ;24 W m 22 is rather larger than estimates from models (e.g., Trenberth et al. 2009;Wild et al. 2015) and the global latent heat flux ;80 W m 22 is lower than some other estimates (e.g., Wild et al. 2013;Stephens et al. 2012). Other fluxes obtained in the NEWS solution are closer to alternative estimates in the literature; for example, the global precipitation is compatible with Adler et al. (2012) and Trenberth et al. (2009), and the global downwelling longwave radiation is similar to those in Stephens et al. (2012), Wild et al. (2013), and Wang and Dickinson (2013). The considerable variation in published values for these global flux components in particular prompts us to explore how these quantities are sensitive to adjustments to the solution method.

Ocean fluxes and their uncertainties
In this section we consider changes to the treatment of the heat fluxes over the oceans in the NEWS solution method. We assess the impact both on the global average flux optimization and on the North Atlantic NSF in particular.

a. Correlated uncertainties in the NEWS solution
We first note that in the original flux data (Fig. 2a) the North Atlantic NSF substantially contrasts with other basins in showing much less heat gain. However, this NSF contrast is not retained in the optimized solution ( Fig. 2b). If all the major basins were adjusted by a more similar amount then the North Atlantic would exhibit a substantial heat loss, more consistent with a strong AMOC transport. To retain the contrast in the North Atlantic fluxes during the optimization we can explicitly allow for spatially correlated uncertainties in the original fluxes, reflecting a modified observation error covariance matrix [R obs in Eq. (6)]. These correlated uncertainties could, for example, be caused by large-scale structural uncertainties in the experimental methods used to derive the regional fluxes from the raw data.
In fact the NEWS solution method is somewhat inconsistent in its approach to spatial error correlations. The monthly, regional, and global uncertainties given in L' Ecuyer et al. (2015, their Table 2) and Rodell et al. (2015a) are only self-consistent if all monthly and regional uncertainties in each flux component are perfectly correlated. However, this would prevent a strong NSF constraint being imposed simultaneously over the multiple land areas. The regional NEWS solution therefore assumes a diagonal (uncorrelated) error covariance matrix when solving in the individual regions, while imposing a global constraint on each flux to match the result obtained using the global annual mean values in L' Ecuyer et al. (2015, their Table 2).
In the following solutions (denoted AllCorr and PLECorr) we aim for a more consistent treatment of uncertainties, listing all of the resulting global fluxes in Table 1. The leftmost result in the table is our implementation of the NEWS solution. Note that in all solutions the TOA flux values are as follows (all in W m 22 ): ISR 5 340.0 6 0.1, OSR 5 99 6 2, and OLR 5 240 6 2. The slightly different OSR compared to the original NEWS value is due to the exclusive use of CERES data in our study.

b. AllCorr solution: All fluxes correlated in major ocean basins
The North Atlantic NSF contrasts shown in Figs. 2a and 2b strongly suggest that spatial error correlations need to be represented explicitly over the oceans to get a more physically reasonable result. Considering Fig. 2a it is clear that the large ocean basins (excepting the North Atlantic) are all initially gaining large amounts of heat and, after optimization, have values of NSF much closer to zero. The small basins (Mediterranean, Arctic Sea, and Black Sea) are not highly out of balance, but the Caribbean has a very large heat gain, which is more similar to the North Pacific. In the AllCorr solution we therefore impose high spatial covariances between prior uncertainties in the North and South Pacific, the North and South Atlantic, the Indian, and the Caribbean basins, for each surface energy flux and precipitation (correlations between different fluxes remain zero). To do this we modify R obs by changing the relevant off-diagonal error covariances such that each has a Pearson correlation coefficient of 0.99, an idealized number that is chosen to impose almost perfect correlation between each region. All other regions are allowed to have spatially uncorrelated errors. Figure 4a shows the optimized NSF for the AllCorr solution. The strongly negative North Atlantic NSF implies that this basin is now losing a large amount of heat to the atmosphere. The values of NSF in the Pacific and Indian basins do not exhibit the large differences observed in the North Atlantic. They are compatible with the values found in the NEWS solution (Fig. 2b) and the ERA-Interim version (Fig. 3). This reflects the fact that only the North Atlantic is an outlier in the input NSF (Fig. 2a) and the AllCorr solution simply preserves its anomalous status using error covariances. The Atlantic Ocean cross-equatorial heat flux implied by the AllCorr solution is 0.52 6 0.14 PW northward, which is much larger than in the NEWS solution. The regional ocean uncertainties are now smaller than those in

c. PLECorr solution: Precipitation and latent heat correlated in major ocean basins
The PLECorr solution is the same as the AllCorr solution except we take solely the latent heat flux and precipitation uncertainties to have high correlations (0.99) between the six larger ocean basins. The surface radiative and sensible heat flux component errors are taken as uncorrelated, enabling a stronger focus on the impact of the correlations on the water cycle. Figure 4b shows the NSF for the PLECorr solution; this solution still allows the North Atlantic to show large net heat losses, although not as large as in Fig. 4a. The Atlantic Ocean cross-equatorial heat flux is equal to 0.34 6 0.43 PW northward. We note that the six correlated basins account for ;85% of the global evaporation, and the larger error covariances for both latent heat and precipitation permit larger global adjustments to the water cycle terms relative to the radiative and sensible heat fluxes. As shown in Table 1 there is now a large increase in global latent heat flux losses. Comparing this result with previous global energy budgets it can be seen that the PLECorr solution is much closer to Stephens et al. (2012), who find strong global latent heat and precipitation fluxes of 88 6 10 W m 22 . As the water cycle adjusts more globally, the other fluxes have adjusted by smaller amounts, including the sensible heat fluxes, which slightly decrease in the PLECorr solution. The regional uncertainties in the PLECorr solution lie between those of the AllCorr solution and the NEWS solution because the implied global uncertainties are only amplified for latent heat and precipitation and not for the other flux terms. The AllCorr and PLECorr solutions both have more realistic-looking ocean NSF than the original NEWS solution and also illustrate the role of spatial error covariances in controlling the global flux solutions. We now consider some evidence for spatial covariability by looking at alternative ocean turbulent flux products. 398 6 2 398 6 2 397 6 2 398 6 2 398 6 2 USW 22 6 1 2 2 6 1 2 2 6 1 2 2 6 1 2 2 6 1 SH 24 6 1 2 4 6 3 2 3 6 1 1 9 6 1 2 4 6 1 LE 80 6 2 8 0 6 3 8 6 6 2 8 5 6 2 8 0 6 2 P 8 0 6 2 8 0 6 3 8 6 6 2 8

d. AltTH solution: Alternative ocean turbulent flux products
The NEWS study used sensible and latent heat fluxes taken from the SeaFlux v1 product (Curry et al. 2004). Given the differences observed in the NEWS latent heat solution relative to other studies, it is interesting to consider whether the use of an alternative flux product would influence the solution and potentially bring it closer to the other estimates.
Tables 2 2017) and Valdivieso et al. (2017). The last row of each table shows the globally averaged fluxes when combined with the land surface fluxes used in the NEWS study. The regional uncertainties on the SeaFlux values are taken from Rodell et al. (2015a), although the global ocean uncertainties assume uncorrelated basins. The mean and standard deviation of the five flux products are also shown in the table. The product standard deviation is generally very consistent with the regional uncertainties provided for SeaFlux, supporting the regional uncertainty estimates used in the NEWS study.
We can also assess interproduct covariability. The mean interbasin correlation coefficient for the six major basins is 0.9 6 0.1 for latent heat and 0.6 6 0.3 for sensible heat fluxes, going some way to justify the assumption of correlated errors between these basins in the AllCorr and PLECorr solutions above.
The J-OFURO v3 product is of particular interest due to its smaller Bowen ratio and larger net turbulent heat losses over the North Atlantic. Compared to SeaFlux v1, the J-OFURO v3 flux quality has been improved, benefitting from its use of multiple satellite sensors with state-of-theart atmospheric humidity retrieval algorithms. In the AltTH solution we replace the SeaFlux v1 turbulent heat fluxes with those from J-OFURO v3, retaining the SeaFlux uncertainties (the interproduct standard deviation could also be used), and without assuming any spatial error covariances. We have, however, retained the SeaFlux Arctic fluxes as the J-OFURO fluxes lack good coverage in that region.

e. Seasonal cycle with correlated ocean fluxes
In this section we consider the seasonal cycle for the PLECorr solution, for which the precipitation and latent heat fluxes are correlated over the six major ocean basins. We retain the monthly SeaFlux uncertainty estimates and use the same spatial error correlations (0.99) each month. The monthly water cycle components for both the PLECorr solution and the NEWS solution, and the difference between the two, are shown in Figs. 6a-c and the monthly NSF is shown in Fig. 6d.
The adjustments represent an intensification of the water cycle; both precipitation ( Fig. 6a) and evaporation (Fig. 6b) increase by ;10% throughout the year in all ocean basins. The changes in atmospheric water convergence ( Fig. 6c) are smaller, with slightly more moisture divergence over the Atlantic basins and more convergence over the Pacific, and almost no impact on land areas. The increased latent heat losses at the ocean surface are partly compensated by other component fluxes in the NSF (Fig. 6d), with all basins showing downward radiation fluxes slightly larger and upward radiation and sensible fluxes slightly smaller. We will compare this seasonal response with that produced by a reanalysis constraint in section 4.

Reanalysis transport constraints
We noted earlier that we had dispensed with the use of the MERRA reanalysis atmospheric water convergences since the surface runoff product supplies the equivalent long-term mean information over land. In this section we investigate reintroducing reanalysis constraints into the inverse model as priors to be adjusted along with the EO-based fluxes.
Adding time-mean horizontal atmospheric energy flux convergences over land would be ineffective because we already impose a strong surface constraint of very small mean NSF in these regions. Similarly, runoff data from land areas offer a sufficient constraint for closing the freshwater budgets and have uncertainties of similar magnitude to those assigned to the MERRA transports in the NEWS study. Over the ocean basins, however, neither (sub)surface energy nor freshwater storage has any regional constraints. Introducing priors for horizontal energy and water convergences from either atmospheric or ocean reanalyses should have an equivalent impact on the regional energy and water budgets.
Here we choose to introduce horizontal heat flux convergences based on ocean reanalysis data rather than use atmospheric reanalysis data for two reasons: first, atmospheric reanalysis energy convergences can be challenging to calculate due to the need to apply a mass correction to the water vapor transport (Mayer et al. 2017); second, atmospheric reanalysis convergences over the oceans have open boundaries at the coasts and will not necessarily sum to small ocean net values (see, e.g., Fig. 2 of Valdivieso et al. 2017). In contrast, ocean reanalysis convergences must closely sum to zero over the global ocean and there is no need to apply a mass correction. For consistency, we use freshwater convergences from the reanalyses that are used to determine the energy convergences; additionally, freshwater convergences from an atmospheric product may conflict with the runoff data provided, as indeed was noted when using MERRA data in Rodell et al. (2015b).  Ferry et al. 2012). Each reanalysis was produced using the ocean model NEMO v3 coupled to the Louvain-la-Neuve sea ice model (LIM) v2 with a model resolution of 0.258. All products assimilated sea surface temperature (in situ observations, satellite data, or both), sea level from satellite altimetry, subsurface temperature and salinity profiles, and sea ice concentration, together with atmospheric forcing taken from ERA-Interim. The heat transport across all sections bounding each of the nine oceanic regions used in the inverse model study is calculated as where r is the density of seawater, c p is the specific heat capacity of seawater, u is the potential temperature, v is the velocity across the section,n is the unit normal vector of the section, and the integral is calculated over the section (Ganachaud and Wunsch 2001). The freshwater transports are calculated as where S is the salinity and S 5 35 psu is the global average ocean salinity (Ganachaud and Wunsch 2001;Wijffels 2001). The heat and freshwater flux convergences in each region are then calculated from each of the four ocean reanalysis products and the mean is determined, with the uncertainty defined as the standard deviation of the products around the mean. When computing the freshwater convergences the runoff data are also accounted for.

a. ORA solution: Transport constraints from ocean reanalyses
In the ORA solution the mean and standard deviation of the heat and freshwater transport constraints are used to include prior values of the variables SEC and SWC [see Eqs.
(2) and (4)] over the oceans. Note this approach is different from the direct combination of CERES TOA fluxes and ERA-Interim convergences (described in section 2c) because the reanalysis constraints are only provided as weak (comparatively uncertain) priors and do not dominate the solution. No error correlations are imposed in this case.
The global fluxes produced in the ORA solution are listed in Table 1. These are unchanged from the NEWS solution as, on a global basis, the convergences of heat and freshwater sum to zero and we have not altered the balance of flux adjustments between the radiative and turbulent contributions at the surface.
The regional NSF incorporating these ocean convergences into the NEWS solution is shown in Fig. 5b. There are several striking features of this solution. The heat loss in the North Atlantic has a much larger central value than in the NEWS fit and although the Caribbean is still gaining heat the magnitude of the increase is greatly reduced from previous solutions. The Arctic now shows surface heat loss, although rather less than in the ERA-Interim constrained solution. This solution corresponds to a cross-equatorial heat transport in the Atlantic Ocean of 0.73 6 0.13 PW northward, very similar to the ERA-Interim derived solution. The ocean heat transport from the North Atlantic into the Arctic is 35 6 5 TW. This is the only solution that produces a realistic value of the Arctic net energy budget (Tsubouchi et al. 2018). Given the sparsity of observations it is not surprising that EO data alone over the Arctic are insufficient to separate the atmospheric and oceanic energy and water transports, even if the TOA data are reliable.
Another feature of this reanalysis-constrained solution is the small regional uncertainties assigned to the optimized ocean NSF values compared to most solutions in section 3. This is the result of introducing a net flux constraint, just as introducing a strong constraint on NSF over land leads to small uncertainties.
The values of NSF in the Pacific and Indian basins computed in the ORA solution are still consistent with those found in the NEWS solution (Fig. 2b) and the ERA-Interim version (Fig. 3).

b. Seasonal cycle results
One way of clearly separating the impact of the EO fluxes and the reanalysis convergences is to solve for the seasonal cycle while only applying the reanalysis convergences as time-mean constraints. Any seasonal variability will thus result purely from the EO flux data. The method is similar to the monthly solution presented in Rodell et al. (2015b), which constrains the time mean fluxes but allows their monthly values to vary. The seasonal results including these annual convergences are shown in Fig. 7, with differences relative to the NEWS solution at the bottom.
The latent heat losses over the North Atlantic increase relative to the NEWS solution throughout the year, but are larger in winter due to the seasonally varying uncertainties. The precipitation also increases over the North Atlantic, reflecting a strengthening of the water cycle over that basin. The other ocean basins show both increases and decreases in latent heat loss, precipitation, and atmospheric water convergence, consistent with no overall global change for the ORA solution in Table 1. The South Pacific exhibits the largest change in atmospheric water convergence, reflecting both increased precipitation and decreased evaporation. In contrast to the PLECorr solution (Fig. 6) the smaller increase in latent heat loss over the North Atlantic, ranging from 3 to 6 W m 22 , is reinforced by changes in the other energy fluxes; the downwelling radiation is reduced and the upward radiation and sensible fluxes are increased, giving a larger net heat loss ranging from 10 to 18 W m 22 and an annual average loss of 16 6 3 W m 22 .

Discussion
The AllCorr solution (section 3a), which retains correlated errors in all fluxes over the six largest ocean basins, might be considered a minimal update to the NEWS solution. The correlations allow the ocean flux adjustments to produce more realistic results while still allowing strong land surface energy constraints to be imposed. Spatial covariances such as those applied in the AllCorr and PLECorr solutions are highly likely to arise in EO-based products due to structural uncertainties in the satellite retrieval methods and instrumental biases. These biases will dominate on larger scales when regional random uncertainties average out (Adler et al. 2012) but it is hard to quantify the spatial error correlations for any single flux product.
We note that simply changing the size of the flux uncertainties without introducing correlations does not impact spatial flux distributions very much. For example, if the default LE and SH uncertainties are doubled (without imposing error correlations) the global and regional fit values are quite similar to those of the NEWS solution albeit with larger uncertainties. The optimized global turbulent heat fluxes in this case are SH 5 27 6 2 W m 22 and LE 5 81 6 2 W m 22 and the North Atlantic net surface flux is 21.4 6 18.7 W m 22 . The lack of correlations means the regional fluxes are not required to adjust in concert and the distinct improvements to the North Atlantic solution are lost.
When considering spatial correlations we have used idealized values of 0.99 but we have also demonstrated that the correlated variability among five turbulent flux products is very high both for the 10-yr means and for the seasonal fluxes in each year (not shown). The interproduct variability does exhibit a seasonal cycle, suggesting that the uncertainties in the different products are dependent on the seasonal environmental conditions. The disadvantage of this method is that the interproduct variability does not account for errors relative to the truth that are shared between the flux products. We have compared the satellite flux products with ICOADS ship data, which can be thought of as a proxy for the truth and, in principle, enables a determination of the true error covariances. The EO product differences relative to ICOADS are correlated on large spatial scales, but the ship data are too spatially and temporally inhomogeneous to take reliable averages for performing a bias correction. Allowing precipitation errors to also be spatially correlated has a significant impact on the solution. If the precipitation errors remain uncorrelated this strongly reduces the adjustments to the water cycle when introducing latent heat error correlations. It is known that the microwave data upon which EO precipitation fluxes are estimated are unable to fully detect weak precipitation (Behrangi et al. 2014) and this bias is likely to be present in all basins and represent a highly correlated error.
We have not sought to investigate spatial error correlations in the radiative flux components although they too could make a significant difference to the results. As noted above it is generally very hard to estimate EO flux errors on large scales because in situ data cannot be used to test error correlations unless they dominate over local temporally varying errors. Improved retrieval methods based on physical processes or state dependence using analyzed atmosphere or ocean conditions are needed. Additionally, correlations between different fluxes (e.g., between precipitation and latent heat flux observations, whose derivations share common data) could also have an important impact. We have tested adding error correlations of 0.99 between the precipitation and latent heat fluxes in each ocean basin in addition to the correlations used in the PLE-Corr solution. The result (not shown) is only slightly different from the PLECorr solution with essentially the same regional and global characteristics. Further data-driven studies to determine additional error correlations may be pursued in the future.
It is also important to note that we have taken the annual and monthly mean regional uncertainties provided by NEWS in Rodell et al. (2015a) as a baseline for our solutions. These uncertainties are only really self-consistent when each component uncertainty is perfectly correlated both spatially and temporally. Therefore by explicitly solving with some uncertainty components and regions correlated and some uncorrelated we reduce the implied global uncertainties relative to those in NEWS. An alternative assumption would be to retain the NEWS regional uncertainties where we intend to explicitly solve with high spatial error correlations, but to inflate regional uncertainties for those flux components to be treated as spatially uncorrelated. This should lead to the same regional distributions of fluxes as we show in the PLECorr solution (e.g., with improved North Atlantic NSF) but with the same global average solution including uncertainties as in the NEWS solution (Table 1).
Of the alternative ocean turbulent flux products we considered, only the J-OFURO v3 product substantially alters both the global average sensible and latent fluxes and the regional net flux solutions over the North Atlantic (the AltTH solution; section 3d). However, even this solution does not correct the North Atlantic fluxes sufficiently to be consistent with the known AMOC heat transports. The J-OFURO v3 product is known to exhibit a dry bias compared to surface observations (Roberts et al. 2019), which partly explains its relatively large evaporation. Other products such as SeaFlux v2 have been developed with the aim of improving fluxes through such an approach (B. Roberts 2018, personal communication) but it is unlikely that all correlated errors can be removed from any of these products over the oceans.
The ORA solution (section 4a), which incorporates ocean reanalysis constraints into the inverse model, is able to effectively provide realistic net surface fluxes over the North Atlantic and the Arctic without having to modify the observation error covariance matrix. When reanalysis transports are used as weak annual mean constraints they do not impact global mean budgets as they only represent a redistribution of time mean surface fluxes, leaving temporal variability driven by EO flux products alone. Furthermore, by removing the use of atmospheric reanalysis transports (MERRA in the original NEWS solutions) we avoid competing with surface runoff information for the annual mean solutions. The seasonal cycle solutions show nonuniform adjustments through the year reflecting seasonal variations in uncertainties. When individual flux component uncertainties are changed, as in the PLECorr solution, the other flux components tend to adjust in the opposite direction, but when net flux constraints are used, as with ocean reanalyses in the ORA solution, all energy fluxes adjust in the same direction. For some basins, in particular the Arctic, exploiting additional ocean data in this way is particularly useful because the narrow oceanic straits greatly reduce uncertainties on inflow/outflow fluxes (Tsubouchi et al. 2018).
Although ocean reanalyses rely on Argo data extending to depths of just 2 km, the ocean heat and freshwater transports are dominated by the upper ocean properties. Small errors in deeper properties will have little impact on these transports. It may be that the surface imbalance needs correcting due to heat changes in the deep ocean but that will not have a big impact on the very large spatial distribution changes shown in this paper.
The Caribbean net surface flux exhibits pronounced variability between the different solutions. This region is too small to receive much correction from the constraint on the global ocean NSF in the uncorrelated solutions. The ocean reanalysis data constrain the Caribbean region to much smaller surface heat gains than the other solutions, demonstrating the value of horizontal (sub)surface transport constraints. More generally, the assumption of a global ocean heat uptake of 1.0 6 0.6 W m 22 does not account for the known variation between different basins. We have tested using the basin-specific values of this uptake that are presented in Cheng et al. (2017). The results change very little because the net basin uptake is always small relative to the flux adjustments and the net basin heat exchanges.
As shown in Table 1, modifying the choice of ocean turbulent heat flux product and/or allowing spatial error covariances to vary for each flux component can have a major impact on the global flux balances while remaining within the range of estimates in the literature. We have adjusted latent and precipitation error covariances together but note that the global water cycle will be more constrained by whichever of these has smaller or less spatially correlated errors. The other global fluxes are very similar to those presented in L' Ecuyer et al. (2015).
We have computed the metric x 2 [ 2J min for all of the fits described in this document. This metric can be used together with the number of degrees of freedom [d, calculated following Ye (1998)] to compute a p value for the hypothesis that the fit values and observations are drawn from the same distribution. The values of x 2 /d for the annual and monthly versions of the PLECorr solution are 14/197 and 236/2679, respectively. Both of these solutions have a p value very close to unity, sustaining the null hypothesis. Similar conclusions hold for all solutions presented in this paper, indicating that the optimized fluxes do not deviate inconsistently from their input values. Indeed, most of the variables move less than one standard deviation from their initial values; the particular exceptions to this occur when a flux has a high error correlation (e.g., precipitation in the PLECorr solution).
With some readjustment of the ocean regions it would be possible to substitute hydrographic observation-based estimates of heat and freshwater transports instead of ocean reanalysis products. The WOCE hydrographic survey estimates (e.g., Ganachaud andWunsch 2000, 2001) could provide example transport products, and time-varying transports might even be used, for example from RAPID (Kanzow et al. 2007;McCarthy et al. 2015) and the Overturning in the Subpolar North Atlantic Program (OSNAP; Lozier et al. 2019) array data. As mentioned above, horizontal transport constraints provide valuable information and will enable the ocean to be divided into smaller regions without making the solutions unconstrained.
The introduction of correlated fluxes and the use of ocean reanalysis constraints have both improved the solution on a global and regional basis. A natural question to ask, then, is what is the optimal way to proceed in future studies. We advocate for using as much information as possible but, critically, the uncertainties associated with EO and reanalysis data must be correctly calculated in order to ensure each component is correctly weighted in the inverse model. If data-driven correlations are included in the observation error covariance matrix, for example, the uncertainty on those correlations should be evaluated and incorporated into the fit.
The time period chosen in this analysis (2000-09) enables a large number of disparate datasets to be employed although data completeness and quality may improve through the period. This can be explored in future studies of interannual variability.

Summary and conclusions
In this paper we have presented several approaches to improve the objective energy and water cycle fluxes first produced in the NEWS study (L'Ecuyer et al. 2015;Rodell et al. 2015b). EO-derived fluxes for TOA radiation and surface energy and water exchanges have been optimally combined under various budget closure constraints for the atmosphere, ocean, and land surface. The input fluxes are adjusted within their uncertainties to produce annual mean and seasonal budgets. Several relative differences to other budgets have been identified in the original NEWS solution; for example, the global sensible (latent) heat fluxes are relatively low (high). There are also strong discrepancies in the net surface energy fluxes over the North Atlantic, Caribbean Sea, and Arctic Ocean that are clearly inconsistent with the fact that the ocean overturning circulation is known to carry substantial heat into these basins. These ocean basins do show much weaker net surface energy gains than other ocean basins in the input EO flux data ( Fig. 2a) but this contrast is lost in the original NEWS optimized solutions (Fig. 2b).
We have shown that substantial improvements can be obtained by imposing error covariances between the larger ocean basins, which share strong NSF imbalances in the original data. As shown in the AllCorr solution (Fig. 4a), if all of the flux components have similarly correlated errors then the regional surface fluxes substantially improve. The NSF in the North Atlantic and Caribbean basins is more consistent with strong AMOC heat transports but the global mean flux components are unchanged. If the flux component error correlations differ then the balance of global flux adjustments alters, as shown in the PLECorr solution (Fig. 4b). We have presented evidence for high error correlations in ocean turbulent heat fluxes and shown that imposing error covariances in the latent and precipitation fluxes alone allows more North Atlantic surface heat losses and a substantial increase in the global net latent heat loss. Replacing the SeaFlux v1 turbulent heat flux product with the more recent J-OFURO v3 also leads to larger global latent heat fluxes, reduced sensible fluxes, and a higher heat loss in the North Atlantic as demonstrated in the AltTH solution (Fig. 5a).
The flux solutions can also be altered by including reanalysis net flux convergences as additional weak constraints within the inverse model (the ORA solution; Fig. 5b). We emphasize this is different from using reanalysis convergences as exact constraints, as imposed by Trenberth et al. (2009), Allan et al. (2014, and Liu et al. (2015). We have chosen to use heat and freshwater convergences over the ocean basins based on four well-studied high-resolution reanalyses, with uncertainties derived from their standard deviation. Over land areas the strong surface flux constraints and runoff data are sufficient to constrain both energy and water cycles. Ocean reanalyses, unlike atmospheric reanalyses, have the advantage of always giving zero global heat flux convergence and a global water convergence made up only of runoff from the land. The ORA solution exhibits improved regional net surface fluxes particularly over the North Atlantic, Caribbean, and Arctic. The balance of the radiative and turbulent heat flux adjustments is not altered, however, and therefore the global latent and sensible fluxes obtained by the original NEWS solution are unchanged. We also derive new seasonal cycle fluxes with ocean reanalysis fluxes only included as an annual mean constraint. We show that the adjustment to the seasonal fluxes is not uniform, being influenced by the seasonal cycle in the original flux uncertainties leading to larger adjustments in wintertime turbulent heat fluxes where the uncertainties are largest.
The approaches outlined in this paper use both EO flux products and their uncertainties along with reanalysis output, with the combination of methods inspired by the vision outlined by the CLIVAR CONCEPT-HEAT program. The results provide adjustments to the NEWS optimized flux solutions as well as showing how different constraints and error assumptions can influence the objective fluxes obtained by the inverse model. There are opportunities to further improve the results by improving the underlying flux products and/or the uncertainties attached to them. The method can also be extended to interannual flux variability, as noted in the original NEWS papers, and further work in this direction should provide valuable solutions for testing the representation of such variability in climate models.

Full Cost Function
We present the most general inverse model solution when some fluxes have prior measurements, some are indirectly derived, and there is a mixture of exact and inexact constraints. Fluxes with (without) prior measurements are denoted F (G) and Lagrange multipliers are contained in the vector l. An example of a flux with (without) a prior measurement is DLR (SEC); an example of an inexact constraint is the global ocean energy balance that is constrained to 1.0 6 0.6 W m 22 [Eq. (2)]; and an example of an exact constraint is the requirement that the atmospheric water cycle must balance to zero in each region [Eq. (3)].
The cost function, introduced in abbreviated form in Eq. (6), is J(F, G, l) 5 1 2 (F 2 F obs ) T R 21 obs (F 2 F obs ) where the first term expresses the offsets between F and the observed values F obs (with error covariance matrix R obs ) and the second and third terms express the K physical balance constraints, both exact and inexact, mentioned in Eq. (6). The quantities v F,G and w F,G are matrices representing linear combinations of the elements of F and G, respectively, V and R V are the central values and covariances of the inexact constraints, and W are the values of the exact constraints. In practice the matrix R V is diagonal so the inexact constraints reduce to the form given in Eq. (7); the term S(F, G) in that expression is equal to v F F 1 v G G.
The derivatives of the cost function with respect to F, G, and l are Given the fact that the cost function is minimized by setting each derivative to zero, these equations can now be re-expressed as follows: Xb 5 a, This equation is solved with least squares techniques. In practice the matrix X is often very ill conditioned and preconditioning must be applied prior to performing the inversion.
The uncertainties on the values of F and G are determined by using the fact that X is equal to the Hessian matrix, = 2 J, and its inverse can be equated with the fit covariance matrix by removing any rows and columns that correspond to Lagrange multipliers. Both correlations in R obs and the presence of exact and inexact constraints influence the structure of the covariance matrix.

Determination of Atlantic Cross-Equatorial Heat Flux
The Atlantic cross-equatorial heat flux can be determined using the fit values of the energy convergence in each ocean basin. Given a vector of ocean SEC values n with covariance matrix R n , the objective is to determine the basin-to-basin transports m. We define a connection matrix C that expresses the links between each basin: it has values of 61 for connected basins, where the sign depends on the direction, and 0 if there is no connection.
The cost function used to determine basin-to-basin transports is J xeq 5 1 2 (n 2 Cm) T R 21 n (n 2 Cm) .
The solution to this equation is simply m 5 C 21 n .
All Southern Hemisphere basins are combined in order to eliminate an ambiguity due to the fact that there are no land barriers in the Antarctic Circumpolar Current. Additionally it is assumed that all transports through the Bering Strait are negligible. These assumptions ensure that the matrix C is full rank and can be inverted without instabilities. The Atlantic cross-equatorial flux is then taken from the relevant element of m. The Hessian matrix of the cost function is equal to CR 21 n C T , which can be inverted to estimate the covariance matrix of the derived transports.