Comparison of ocean heat content from two eddy-resolving hindcast simulations with OFES1 and OFES2

The ocean heat content (OHC) estimates from high-resolution hindcast simulations from the Ocean General 10 Circulation Model for the Earth Simulator Version 1 (OFES1) and Version 2 (OFES2), and a global objective analysis of subsurface temperature observations (EN4.2.1) were compared. There was an OHC increase in most of the global ocean over a 57-year period, mainly a result of vertical displacements of neutral density surfaces. However, we found substantial differences in the temporal and meridional distributions of the OHC between the two OFES hindcasts. The spatial distributions of potential-temperature change also differed significantly, especially in the Atlantic Ocean. The spatial distributions of the 15 time-averaged surface heat flux and heat transport from the OFES1 and OFES2 were highly correlated, but differences could be seen. However, these differences, more specifically in the heat transport, were only partially responsible for the OHC differences. The marked OHC differences may arise from the different vertical mixing schemes and may impact the largescale pressure field, and thus the geostrophic current. The work here should be a useful reference for future OFES users.

dissipation (Jayne and St. Laurent 2001;St. Laurent et al., 2002), whereas the OFES1 uses the K-profile parameterization scheme (Large et al., 1994). With the oceanic field on 1 st Jan 1958 from the OFES1 as the initial conditions, the OFES2 used here has been integrated from 1958 to 2016. To limit the computation cost, we subsampled the OFES1 and OFES2 simulations every 5 grid points in the horizontal direction. 85 To validate the OHC from the two sets of OFES data, we used the EN4 from the UK Met Office Hadley Centre as a reference. The monthly EN4 data can be considered as an objective analysis that is primarily based on observations (Good et al., 2013), with a horizontal resolution of 1° and 42 vertical levels from 5 m to 5350 m. The EN4 assimilates data mainly from the World Ocean Database (WOD) and the Coriolis dataset for ReAnalysis (CORA). Pre-processing and quality checks are conducted before the observational data are used to construct this objective analysis product. 90 Although we use the EN4 results as a reference for evaluating the OFES performance in simulating the 57-year ocean thermal state, it is also worthy to note that the EN4 cannot be taken as the actual ocean state. The main reason is that the measurements used to construct the EN4 datasets are sparse and inhomogeneous in both the temporal and spatial domains, and far from sufficient to resolve mesoscale or even sub-mesoscale motions. There were more observations in the northern hemisphere than in the southern hemisphere, and a higher density of records became available only after the World Ocean 95 Circulation Experiment (WOCE) in the 1990s and installation of the Argo profiling floats in the 2000s. Table 1 summarizes these three ocean datasets.

Methods
We compared the three datasets during the period between 1960-2016. Although the two OFES datasets may not be well spun 100 up in the beginning, especially the OFES2, we focus on their differences on a multi-decadal scale. Following convention, the OHC values here are the OHC anomalies relative to 1960. At each grid point, the OHC was calculated as where is the seawater density, the grid volume, the specific heat of seawater at constant pressure, the potential temperature and 1960 the potential temperature in 1960. A value of 4.1×10 6 kg J m -3 K -1 was used for the product of and 105 specific heat of seawater (Palmer et al., 2011).
Both the global and individual-basin OHCs were calculated for comparison. Figure 1 shows the domains of the Pacific, Atlantic and Indian Oceans between 64°S and 64°N, with their respective marginal seas included.  The changes in the potential temperature Δθ were decomposed into an SP component (second term in Eq. (2) below) and an HV component (third term in Eq. (2)) (Bindoff and McDougall 1994). HV is the Eulerian measure of Δθ at fixed depths, resulting from the vertical displacement of neutral surfaces (Häkkinen et al., 2016). SP represents changes along the neutral surfaces. This decomposition helps to identify the dominant mechanisms changing the potential temperature. The formula 115 decomposing the potential temperature is We used the program by Jackett and McDougall (1997) to calculate the neutral densities, HV and SP. The main inputs for this program are the potential temperature θ and salinity S. As the code limits the latitude domain to between 80°S and 64°N, we set our investigation domain to between 64°S and 64°N to avoid comparisons in sea-ice impacted areas (only OFES2 120 includes a sea-ice model).

Results
The principal aim here is to compare the results from the OFES1and OFES2, with the EN4 acting as an observation-based reference. If there is a significant difference between the OFES2 result and that of one or both of the other two datasets, does this represent a real phenomenon not present in the other two widely used datasets or is it an unwanted property of the newly 125 released OFES2 simulation? In this section, we compare the three sets of results for the global ocean, and for each of the Pacific, Atlantic and Indian Oceans individually. The first task is to identify significant differences.

Time evolution of the OHC, HV and SP from 1960 to 2016
In this section, we compare the OHC time series from the three datasets and with some important findings in the literature. In Fig. 2a, the EN4 shows that the global upper ocean experienced cooling during the 1960s, followed by an approximately linear warming ("linear warming" used here is as a short-hand for "warming at a linear rate") since the 1970s (Achutarao et al., 2007;Zanna et al., 2019). Although this cooling is reproduced in both the OFES datasets, the linear warming appeared in rate in the global upper-layer ocean after 1994 was around 4.3 ZJ/yr from the EN4 and 4.9 ZJ/yr from the OFES1 (1 ZJ = 10 21 J 135 and yr means year), but only 1.0 ZJ/yr from the OFES2. Unlike in the other datasets, a sharp and remarkable OHC reduction stood out in the OFES1 in 1987 (Fig. 2a).
In the upper layer of the Pacific Ocean (Fig. 2b), The result from the EN4 shows that the upper layer of the Pacific Ocean was largely warming with some sporadic exceptions. Conversely, the OFES products indicate an overall cooling but in different ways. More specifically, the OFES1 indicates cooling before 1987 in the upper Pacific Ocean, with a sudden cooling 140 https://doi.org/10.5194/gmd-2021-95 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License. then occurring. The cooling trend then reversed to warming. In the upper layer of the Atlantic Ocean (Fig. 2c), the OHC time series from the EN4 and the OFES1 are highly correlated (Fig. 2c), with a correlation coefficient of 0.9. In addition, the overall warming rate was around 1.0 ZJ/yr from the EN4 and 0.8 ZJ/yr from the OFES1. Strikingly, the OFES2 presents a notable cooling of 1.4 ZJ/yr before 2000. Overall, the absolute differences between the three products in the upper layer of the Indian Ocean are the smallest (Fig. 2d), with the OFES2 showing larger OHC increase before 2010. A summary of the total warming 145 rate in the upper ocean from 1960-2016 is given in Tab. 2. When compared to the EN4, the OFES data show much weaker warming or even cooling, especially the OFES2.
The calculated HV and SP can help identify how the ocean water warms or cools (Häkkinen et al., 2016). The HV dominated the OHC variations (Figs. 2e-h) and evolved substantially different to the SP (Figs. 2i-l). This means that the OHC variations largely result from the vertical movement of the neutral surfaces over this 57-year period, as seen in all three products. The 150 HV dominance in OHC variations was also examined in Häkkinen et al. (2016). One interesting point is that the EN4 and

Meridional distribution of the zonal-integrated OHC, HV and SP
The resolutions of the EN4 and the OFES hindcasts are different, we calculated the OHC increase in each latitude band of 1° from 64°S to 64°N, rather than on the native grid of these data. For clarity of presentation, we denoted these two periods as P1 (1960-1962) and P2 (2014-2016).

a The upper ocean (0-500 m)
In the upper layer of the global ocean (Fig. 4a), both the EN4 and OFES1 show an almost consistent warming in the upper 200 extratropical region, except for the region 60° poleward (Fig. 4a). Moreover, these two datasets capture a common warming peak of approximately 7 ZJ at around 40°S. The major differences between the OFES1 and the EN4 are seen in the upper tropical ocean, a region with consistent warming in the EN4 but two troughs at around 15°S and 15°N in the OFES1. The two OFES datasets give similar trends in the southern tropical region, with a trough of OHC indicated by both.
In each of the major basins, the EN4 shows moderate OHC increases in the upper layer for almost the whole examined 205 latitude range, with a few sporadic exceptions. In the upper layer of the Pacific Ocean (Fig. 4b), the two OFES datasets are more similar to each other in the Southern Hemisphere than in the Northern Hemisphere. Similar to the EN4, the OFES1 also indicates that warming existed over almost the whole latitude range of the upper layer of the Atlantic Ocean (Fig. 4c). In contrast, the OFES2 suggests marked OHC decreases between 50°S and 30°N of the upper Atlantic Ocean. In the upper layer of the Indian Ocean (Fig. 4d), all three datasets present quite small OHC variations in the Northern Hemisphere. However, 210 south of the Equator, the OFES1 is highly close to the EN4, both of which presenting a notable warming peak at around 45°S.
A distinct warming peak also appears in the OFES2, but about 15° equatorward. The HV clearly contributes the bulk of the  In the global intermediate ocean (Fig. 5a), the meridional OHC distributions show again the closeness between the EN4 and 220 OFES1, both having a marked warming peak between around 40-45°S. Their major disagreements are seen in the northern tropic. Between 40°S and 30°N, there are negative OHC variations in the OFES2, which also missed the large warming peak (between around 40-45°S) above-mentioned by the other two.
As for the meridional OHC distribution in the intermediate layer of the Pacific Ocean (Fig. 5b), the OFES2 agrees well with the EN4 between around 64°S and 25°S. To the north of 60°S in the intermediate layer of the Pacific Ocean, the OFES1 225 indicates a consistent OHC increase, whereas the OHC is shown to slightly decrease between 30°S and 20°N in the OFES2.
As indicated by the EN4, warming is present at almost all the latitudes considered in the intermediate layer of the Atlantic Ocean (Fig. 5c). However, cooling can be seen in the OFES data in different latitudal intervals. The OFES1 resembles the EN4 in the southern Indian Ocean (Fig. 5d) in presenting a significant warming peak slightly to the north of 50°S, but this peak is missed by the OFES2. In spite of this distinction, the OFES2 also reveals a large OHC increase near the Equator, just as the 230 other two.

Spatial patterns of the potential-temperature changes
The two previous sections described the global and basin-wide OHC distributions in the temporal domain and the longitudinal direction quantitatively. To further investigate the detailed agreements and discrepancies in the warming or cooling from these three datasets, we calculated the volume-averaged potential temperature θOHC. This was calculated by dividing the total OHC 240 variations in the water column of each ocean layer (upper or intermediate) by the corresponding total water volume and the product of seawater density and specific heat capacity (ρ×Cp = 4.1×106 kg Jm -3 K -1 ). We then compared the change in the volume-averaged potential temperature, ΔθOHC (P2 average (2014-2016) minus P1 average (1960)(1961)(1962)) from the different datasets, as shown in Fig. 6. We also calculated ΔθHV and ΔθSP, derived from the HV and SP, respectively, in a similar way.
The reason for using Δθ rather than the OHC is that the latter is an extensive quantity, but its variation at each grid is directly 245 related to the Δθ. To facilitate interpreting the results, we defined major water masses for both the upper and intermediate oceans in each basin, following Emery (2001), as shown in Tabs. 4 and 5. Readers are referred to Emery (2001) for more details. The geographical distribution of the major water masses analysed here can be found in Emery (2001). In the intermediate layer of the Pacific Ocean, a clear difference between the OFES1 and others is that the former presented most of the Pacific Subarctic Intermediate Water (PSIW) was warming (Fig. 7b). Both the EN4 and OFES1 indicate a slight warming in the California Intermediate Water (CIW), but cooling is shown for the CIW in the OFES2. All three datasets consistently suggest warming in the southern flank of the Antarctic Intermediate Water (AAIW) in the Pacific Ocean (Figs.  305 7a-c), with the warming shown by the OFES1 strongest. Substantial discrepancies exist in the central and northern parts of the AAIW. More specifically, there is a vast region of cooling in the central part of the AAIW in the southern Pacific Ocean from both the EN4 and OFES2 (Fig. 7a,c), whereas the OFES1 reveals a largely warming pattern. The OFES2 also showed negative Δθ in the region west of the northwestern African continent (Fig. 7c), which cannot be seen in the other two data.
The EN4 indicates that almost the whole intermediate layer of the Atlantic Ocean was warming (Fig. 7a) but there exist 310 regions of intense cooling in the OFES data. For instance, both the OFES1 and OFES2 display cooling in the bulk of the Western Atlantic Subarctic Intermediate Water (WASIW), especially the OFES2 (Fig. 7c). There is a similar warming tongue between 30-60°N in the middle of the Atlantic Ocean in the OFES1 and OFES2, but that in the latter much stronger (Figs. 7b,c). Both the OFES1 and OFES2 show that the southern flank of the Atlantic branch of the AAIW was warming, but the warming in the OFES1 extended over a longer distance northward to around 30°S (Fig. 7b). On the contrary, the EN4 trend 315 on the southern side of the AAIW is much weaker compared to others.
In the Indian Ocean, large discrepancies appear in the Red Sea-Persian Gulf Intermediate Water (RSPGIW): the EN4 shows weak warming (Fig. 7a), whereas the OFES1 and OFES2 indicate intense cooling (moderate warming) and intense warming (moderate cooling) for the eastern (western) parts of the RSPGIW, respectively (Figs. 7b,c). Large differences also happen to in the Indonesian Intermediate Water (IIW). Weak and strong warming is shown by the EN4 and OFES1 (Figs. 7a,b), 320 respectively, but significant cooling in the OFES2 (Fig. 7c). The EN4 and OFES1 display strong warming in the ACC path, whereas intense cooling is indicated by the OFES2. At around 60°S in the Indian Ocean, both the OFES datasets show strong cooling (Figs. 7b,c), but ∆θ is quite small in the EN4.

b The intermediate ocean (500-1400 m)
From the T-S diagram (Fig. 9), it can be seen that there are significant differences in salinity in the OFES1 compared to the EN4 and OFES2, the latter two being similar (for example, Figs. 9a,c). Most of the waters have a warming trend in the EN4 355 and OFES1 (Figs. 9a,b). However, a remarkable OHC decrease can be seen in the OFES2 (Fig. 9c). More specifically for the OFES2, the AAIW is responsible for the bulk of the OHC increase in the Pacific Ocean (Fig. 9g). In the Atlantic Ocean, the OFES2 shows that almost all the OHC increase came from the EASIW and MW (Fig. 9k). The major OHC increase in the Indian Ocean can be attributed to the RSPGIW in the OFES2 (Fig. 9o). The AAIW of the Indian Ocean is found to account for large OHC increase in the EN4 and OFES1 (Fig. 9 m,n), but just opposite in the OFES2 (Fig. 9o).

Analysis of surface heat flux and heat transport
The above analyses showed the differences in the estimated OHC and its related potential temperature change from the OFES1 and OFES2 by comparing them with the observation-based EN4; this section analyses possible causes for the differences. We 365 firstly compare the time-averaged heat flux through the ocean surface from the OFES1 and OFES2 (the source term in the governing equation of temperature; not applicable to the EN4). We also calculate and compare the horizontal and vertical heat transports (the advection terms in the governing equation of temperature), whose differences mainly result from the simulated circulations and temperature fields.
The surface heat-flux patterns from the OFES1 and OFES2 are spatially similar (Figs. 10a, b), but certain differences can 370 be seen in Figs. 10 c,d. Most of the Pacific and Indian Oceans gained less heat or lost more heat in the OFES2, a major exception being the eastern Pacific Equator. Comparing the heat flux distributions in Fig. 10 with the Δθ distributions from the OFES1 and OFES2 in Fig. 6, there is no clearly positive correlation between the surface heat flux and the Δθ. For instance, the more intense heat flux in the OFES2 does not lead to a stronger OHC increase in the eastern Pacific Equator. This implies that the examined OHC or potential temperature differences are caused by other factors. 375 https://doi.org/10.5194/gmd-2021-95 Preprint. Discussion started: 20 April 2021 c Author(s) 2021. CC BY 4.0 License.
The large-scale distribution of the horizontal heat transport is a direct reflection of the general circulation and the upper circulation is mainly wind driven. For simplicity, we will not describe the general patterns of heat transport, which is highly similar from these two data, as can be seen from Fig. 11. Instead, we will mainly focus on the differences between the OFES1 395 and OFES2 aiming at explaining the OHC and ∆θ disagreements shown above. Specifically, the stronger westward heat transport in the vicinity of the Pacific Equator (Fig. 11c) may be responsible for the weaker cooling associated with the southern flank of the WNPCW and weaker warming of the eastern Pacific Equator (Fig. 6c). The larger eastward heat transport associated with the Kuroshio can also partially explain the weaker warming of the WNPCW within the Kuroshio pathway in the OFES2, at least locally. The stronger ZHT and MHT associated with the Gulf Stream (Figs. 11c, f) in the OFES2 lead to 400 more heat being transported out, and thus to the local weak cooling in Fig. 6c, compared to moderate warming in the OFES1 (Fig. 6b). The notable warming tongue in the middle of the upper Atlantic Ocean at around 35-64°N may result from the stronger MHT in the OFES2 (Fig. 11f). However, the extensive cooling trend in the vast area of the upper Atlantic Ocean cannot be well explained by the heat advection here and heat flux in Fig. 10, and thus may be attributed to heat diffusion, which cannot be examined here. In the upper layer of the Indian Ocean, the larger amount of heat transported westward in the 405 southern tropics may be responsible for the stronger cooling trend in the IUW and the much weaker cooling in the western flank of the IEW in the OFES2. In addition, the ZHT associated with the northern part of the ACC in the Indian Ocean is stronger in the OFES2, which may account for the cooling there. The much stronger warming of the central south Indian Ocean may be a result of the lesser downward heat transport there in the OFES2 (Fig. 11i).