Ocean eddies play an important role in the distribution of heat, salt, and other tracers in the global ocean. But while surface eddies have been studied extensively, deeper eddies are less well understood. Here we study deep coherent vortices (DCVs) in the Northeast Atlantic Ocean using a high resolution numerical simulation. We perform a census of the DCVs on the $27.60$ kg/m$^3$ isopycnal, at the depth of $700-1500$ m, where DCVs of Mediterranean water (meddies) propagate. We detect a large number of DCVs, with maxima around continental shelves, and islands, dominated by small and short-lived cyclones. However, the large and long-lived DCVs are mostly anticyclonic. Among the long-lived DCVs, anticyclonic meddies, stand out. They grow in size by merging with other anticyclonic meddies. Cyclonic meddies are also regularly formed, but most of them are destroyed near their formation sites due to the presence of the energetic anticyclonic meddies, which destroy cyclones by straining and wrapping the positive vorticity around their core. During their life cycle, as they propagate to the southwest, anticyclonic meddies can interact with other DCVs, including anticyclones containing Antarctic Intermediate Water generated near the Moroccan coast, Canary anticyclonic DCVs and cyclonic DCVs generated south of $30^\circ$N along the African continental shelf. With these latter, they can form dipoles, and with the former, they co-rotate pro tempore. Thus, a more detailed view of the life cycle of anticyclonic meddies is proposed: they grow by merging, undergo multiple interactions along their path, and they decay at low latitudes.

We use seismic ambient noise recorded by dense ocean bottom nodes (OBNs) in the Gorgon gas field, Western Australia, to compute time-lapse seafloor models of shear-wave velocity. The extracted hourly cross-correlation (CC) functions in the frequency band 0.1 – 1 Hz contain mainly Scholte waves with very high signal to noise ratio. We observe temporal velocity variations (dv/v) at the order of 0.1% with a peak velocity change of 0.8% averaged from all station pairs, from the conventional time-lapse analysis with the assumption of a spatially homogeneous dv/v. With a high-resolution reference (baseline) model from full waveform inversion of Scholte waves, we present an elastic wave equation based double-difference inversion (EW-DD) method, using arrival time differences between the reference and time-lapsed Scholte waves, for mapping temporally varying dv/v in the heterogeneous subsurface. The time-lapse velocity models reveal increasing/decreasing patterns of shear-wave velocity in agreement with those from the conventional analysis. The velocity variation exhibits a ~24-hour cycling pattern, which appears to be inversely correlated with sea level height, possibly associated with dilatant effects for porous, low-velocity shallow seafloor and rising pore pressure with higher sea level. This study demonstrates the feasibility of using dense passive seismic surveys for quantitative monitoring of subsurface property changes in the horizontal and depth domain.

This manuscript reports on a robot called Clio that we developed to facilitate basin-scale studies of ocean microbial communities and their biochemistry, to better understand how marine microorganisms regulate ocean and Earth system environmental cycles. Clio is designed to facilitate global-scale studies of ocean biochemistry, to move vertically through the water column with high precision, and specifically to return sensor data and samples from large swaths of the ocean ranging in depths from the surface to 6,000 m. Clio is capable of flexible, precise vertical motion that few other ocean robots can perform, and none to our knowledge over this depth range. We tested Clio extensively over several years, six cruises, and 26 dives, it is now fully operational and this manuscript describes all that we did to convince ourselves this was so. In June 2019, it completed its first large-scale ocean survey, and for which this manuscript will be the first data presentation.

Convergent coastal-plain estuaries have been shortened by dam-like structures worldwide. We used 31 long-term water level stations and a semi-analytical tide model to investigate the influence of a dam and landward-funneling on tides and storm surge propagation in the greater Charleston Harbor region, South Carolina, where three rivers meet: the Ashley, Cooper, and Wando. Our analysis shows that the principle tidal harmonic (M2), storm surge, and long-period setup-setdown (~4–10 days) propagate as long waves with the greatest amplification and celerity observed in the M2 wave. All waves attenuate in landward regions, but, as they approach the dam on the Cooper River, a frequency dependent response in amplitude and phase progression occurs. Dam-induced amplification scales with wave frequency, causing the greatest amplification in M2 overtides. Model results show that funneling and the presence of a dam amplify tidal waves through partial and full reflection, respectively. The different phase progression of these reflected waves, however, can ultimately reduce the total wave amplification. We use a friction-convergence parameter space to demonstrate how amplification is largest for partial reflection, when funneling and wave periods are not extreme (often the case of dominant tides), and for full reflection, when funneling and/or wave periods are small. The analysis also shows that in the case of long period events (>day), such as storm surges, dams may attenuate the wave in funneling estuaries. However, dams may amplify the most intense storm surges (short, high) more than funneling with unexpected consequence that can greatly increase flood exposure.

North Atlantic sea surface temperatures (NASST), particularly in the subpolar region, are among the most predictable locations in the world’s oceans. However, the relative importance of atmospheric and oceanic controls on their variability at multidecadal timescales remain uncertain. Neural networks (NNs) are trained to examine the relative importance of oceanic and atmospheric predictors in predicting the NASST state in the Community Earth System Model 1 (CESM1). In the presence of external forcings, oceanic predictors outperform atmospheric predictors, persistence, and random chance baselines out to 25-year leadtimes. Layer-wise relevance propagation is used to unveil the sources of predictability, and reveal that NNs consistently rely upon the Gulf Stream-North Atlantic Current region for accurate predictions. Additionally, CESM1-trained NNs do not need additional transfer learning to successfully predict the phasing of multidecadal variability in an observational dataset, suggesting consistency in physical processes driving NASST variability between CESM1 and observations.

A document by Roy Barkan. Click on the document to view its contents.

Sea ice-waves interactions have been widely studied in the marginal ice zone, at relatively low wind speeds and wave frequencies. Here, we focus on very different conditions typical of coastal polynyas: extremely high wind speeds and locally-generated, short, steep waves. We overview available parameterizations of relevant physical processes (nonlinear wave-wave interactions, energy input by wind, whitecapping and ice-related dissipation) and discuss modifications necessary to adjust them to polynya conditions. We use satellite-derived data and spectral modelling to analyze waves in ten polynya events in the Terra Nova Bay, Antarctica. We estimate the wind-input reduction factor over ice in the wave-energy balance equation at 0.56. By calibrating the model to satellite observations we show that exact treatment of quadruplet wave-wave interactions (as opposed to the default Discrete Interaction Approximation) is necessary to fit the model to data, and that the power n>4 in the sea-ice source term S_ice~f^n (where f denotes wave frequency) is required to reproduce the observed very strong attenuation in spectral tail in frazil streaks. We use a very-high resolution satellite image of a fragment of one of the polynyas to determine whitecap fraction. We show that there are more than twofold differences in whitecap fraction over ice-free and ice-covered regions, and that the model produces realistic whitecap fractions without any tuning of the whitecapping source term. Finally, we estimate the polynya-area-integrated wind input, energy dissipation due to whitecapping, and whitecap fraction to be on average below 25%, 10% and 30%, respectively, of the corresponding open-water values.

Mesoscale eddies play an important role in transporting water properties, enhancing air-sea interactions, and promoting large-scale mixing of the ocean. They are generally referred to as “coherent” structures because they are organized, rotating fluid elements that propagate within the ocean and have long lifetimes (months or even years). Eddies have been sampled by sparse in-situ vertical profiles, but because in-situ ocean observations are limited, they have been characterized primarily from satellite observations, numerical simulations, or relatively idealized geophysical fluid dynamics methods. However, each of these approaches has its limitations. Many questions about the general structure and “coherence” of ocean eddies remain unanswered. In this study, we investigate the properties of 7 mesoscale eddies sampled with relative accuracy during 4 different field experiments in the Atlantic. Our results suggest that the Ertel Potential Vorticity (EPV) is a suitable parameter to isolate and characterize the eddy cores and their boundaries. The latter appear as regions of finite horizontal extent, characterized by a local extremum of the vertical and horizontal components of the EPV. These are found to be closely related to the presence of a different water mass in the core (relative to the background) and the steepening of the isopycnals due to eddy occurrence and dynamics. Based on these results, we propose a new criterion for defining eddies. We test our approach using a theoretical framework and explore the possible magnitude of this new criterion, including its upper bound.

A significant shift in Earth’s climate characterizes the Neogene, transitioning from a single-ice-sheet planet to the current bipolar configuration. This climate evolution is closely linked to changing ocean currents, but globally-distributed continuous high-resolution sedimentary records are needed to fully capture this interaction. The Ocean Drilling Program (ODP) Site 752, located on Broken Ridge in the Indian Ocean, provides such a Miocene-to-recent archive. We use X-ray fluorescence (XRF) core scanning to build an eccentricity-tuned age-depth model and reconstruct paleoceanographic changes since 23 Ma. We find two intervals of enhanced productivity, during the early and middle Miocene (18.5 – 13.7 Ma) and late Pliocene/early Pleistocene (3 – 1 Ma). We also report a mixed eccentricity-obliquity imprint in the XRF-derived paleoproductivity proxy. In terms of grain size, three coarsening steps occur between 19.2 – 16 Ma, 10.8 – 8 Ma, and since 2.6 Ma. The steps respectively indicate stronger current winnowing in response to vigorous Antarctic Intermediate Water flow over Broken Ridge in the early Miocene, the first transient onset of Tasman Leakage in the Late Miocene, and the intensification of global oceanic circulation at the Plio-Pleistocene transition. High-resolution iron and manganese series provide a detailed Neogene dust record. This study utilized a single hole from an ODP legacy-site. Nevertheless, we managed to provide novel perspectives on past Indian Ocean responses to astronomical forcing. We conclude that Neogene sediments from Broken Ridge harbor the potential for even more comprehensive reconstructions. Realizing this potential necessitates re-drilling of these sedimentary archives utilizing modern drilling strategies.

An ensemble of experiments based on a ¼° global NEMO configuration is presented, including tidally forced and non-tidal simulations, and using both the default z* geopotential vertical coordinate and the z~ filtered Arbitrary Lagrangian-Eulerian coordinate, the latter being known to reduce numerical mixing. This is used to investigate the sensitivity of numerical mixing, and the resulting model drifts and biases, to both tidal forcing and the choice of vertical coordinate. The model is found to simulate an acceptably realistic external tide, and the first-mode internal tide has a spatial distribution consistent with estimates from observations and high-resolution tidal models, with vertical velocities in the internal tide of over 50 meters per day. Tidal forcing with the z* coordinate increases numerical mixing in the upper ocean between 30°S and 30°N where strong internal tides occur, while the z~ coordinate substantially reduces numerical mixing and biases in tidal simulations to levels below those in the z* non-tidal control. The implications for the next generation of climate models are discussed.

We explore the sensitivity of Southern Ocean surface and deep ocean temperature and salinity biases in the FOCI coupled climate model to atmosphere-ocean coupling time step and to lateral diffusion in the ocean with the goal to reduce biases common to climate models. The reference simulation suffers from a warm bias at the sea surface which also extends down to the seafloor in the Southern Ocean and is accompanied by a too fresh surface, in particular along the Antarctic coast. Reducing the atmosphere-ocean coupling time step from 3 hours to 1 hour results in increased sea-ice production on the shelf and enhanced melting to the north which reduces the fresh bias of the shelf water while also strengthening the meridional density gradient favouring a stronger Antarctic Circumpolar Current (ACC). With the shorter coupling step we also find a stronger meridional overturning circulation with more upwelling and downwelling south and north of the ACC respectively, as well as a reduced warm bias at almost all depths. Tuning the lateral ocean mixing has only a small effect on the model biases, which contradicts previous studies using a similar model configuration. We note that the latitude of the surface westerly wind maximum has a northward bias in the reference simulation and that this bias is unchanged as the surface temperature and sea-ice biases are reduced in the coupled simulations. Hence, the surface wind biases over the Southern Hemisphere midlatitudes appear to be unrelated to biases in sea-surface conditions.

Recent marine heatwaves in the Gulf of Alaska have had devastating and lasting impacts on species from various trophic levels. As a result of climate change, total heat exposure in the upper ocean has become longer, more intense, more frequent, and more likely to happen at the same time as other environmental extremes. The combination of multiple environmental extremes can exacerbate the response of sensitive marine organisms. Our hindcast simulation provides the first indication that more than 20 % of the bottom water of the Gulf of Alaska continental shelf was exposed to quadruple heat, positive [H+], negative Ωarag, and negative [O2] compound extreme events during the 2018-2020 marine heat wave. Natural intrusion of deep and acidified water combined with the marine heat wave triggered the first occurrence of these events in 2019. During the 2013-2016 marine heat wave, surface waters were already exposed to widespread marine heat and positive [H+] compound extreme events due to the temperature effect on the [H+]. We introduce a new Gulf of Alaska Downwelling Index (GOADI) with short-term predictive skill, which can serve as indicator of past and near-future positive [H+], negative Ωarag, and negative [O2] compound extreme events on the shelf. Our results suggest that the marine heat waves may have not been the sole environmental stressor that led to the observed ecosystem impacts and warrant a closer look at existing in situ inorganic carbon and other environmental data in combination with biological observations and model output.

In the northern Gulf of Mexico shelf, the Mississippi-Atchafalaya River System (MARS) impacts the carbonate system by delivering freshwater with a distinct seasonal pattern in both total alkalinity (Alk) and dissolved inorganic carbon (DIC), and promoting biologically-driven changes in DIC through nutrient inputs. However, how and to what degree these processes modulate the interannual variability in calcium carbonate solubility have been poorly documented. Here, we use an ocean-biogeochemical model to investigate the impact of MARS’s discharge and chemistry on interannual anomalies of aragonite saturation state (ΩAr). Based on model results, we show that the enhanced mixing of riverine waters with a low buffer capacity (low Alk-to-DIC ratio) during high-discharge winters promotes a significant ΩAr decline over the inner-shelf. We also show that increased nutrient runoff and vertical stratification during high-discharge summers promotes strong negative anomalies in bottom ΩAr, and less intense but significant positive anomalies in surface ΩAr. Therefore, increased MARS discharge promotes an increased frequency of suboptimal ΩAr levels for nearshore coastal calcifying species. Additional sensitivity experiments further show that reductions in the Alk-to-DIC ratio and nitrate concentration from the MARS significantly modify the simulated ΩAr spatial patterns, weakening the positive surface ΩAr anomalies during high-discharge summers or even producing negative surface ΩAr anomalies. Our findings suggest that riverine water carbonate chemistry is a main driver of interannual variability in ΩAr over river dominated ocean margins.

Mixing along isopycnals plays an important role in the transport and uptake of oceanic tracers. Isopycnal mixing is commonly quantified by a tracer diffusivity. Previous studies have estimated the tracer diffusivity using the rate of dispersion of surface drifters, subsurface floats, or numerical particles advected by satellite-derived velocity fields. This study shows that the diffusivity can be more efficiently estimated from the dispersion of coherent mesoscale eddies. Coherent eddies are identified and tracked as the persistent sea surface height extrema in both a two-layer quasigeostrophic (QG) model and an idealized primitive equation (PE) model. The Lagrangian diffusivity is estimated using the tracks of these coherent eddies and compared to the diagnosed Eulerian diffusivity. It is found that the meridional coherent eddy diffusivity approaches a stable value within about 20–40 days in both models. In the QG model, the coherent eddy diffusivity is a good approximation to the upper-layer tracer diffusivity in a broad range of flow regimes, except for small values of bottom friction or planetary vorticity gradient, where long-range correlations between same-sign eddies become important. In the PE model, the tracer diffusivity has a complicated vertical structure and the coherent eddy diffusivity is correlated with the tracer diffusivity at the e-folding depth of the energy-containing eddies where the intrinsic speed of the coherent eddies matches the rms eddy velocity. These results suggest that the oceanic tracer diffusivity at depth can be estimated from the movements of coherent mesoscale eddies, which are routinely tracked from satellite observations.

The effects of contemporary increases in riverine freshwater into the Arctic Ocean are estimated from ocean model simulations, using two runoff data sets. One runoff data set which is based on older climatological data, which has no inter-annual variability after 2007 and as such does not represent the observed increases in river runoff into the Arctic. The other data set comes from a hydrological model developed for the Arctic drainage basin, which includes contemporary changes in the climate. In the pan-Arctic this new data set represents an approximately 11% increase in runoff, compared with the older climatological data. Comparing two ocean model runs forced with the different runoff data sets, overall changes in different freshwater markers across the basin were found to be between 5-10%, depending on the area investigated. The strongest increases were seen from the Siberian rivers, which in turn caused the strongest freshening in the Eastern Arctic.

The Beaufort Sea High is a high-pressure system located in the Beaufort Sea and influences ocean circulation in the western Arctic known as the Beaufort Gyre. Wrangel Island, located in the western Chukchi Sea, typically experiences easterly sea ice motion due to the Beaufort Gyre. We find that under these climatological conditions, moving ice is blocked by the island and piles up on its eastern side, while ice on its western side continues to drift. This results in an ice thickness dipole across the island. A reversal in the sense of the oceanic and atmospheric circulation across the western Arctic results in a dipole with the opposing sign. We find the dipole is present throughout the year and is strongest in January when the ice thickness difference is approximately 1m. During the spring, it is associated with the transient opening of a polynya to the west of the island. The dipole is the result of opposing ice divergence and convergence across the western Arctic and may impact ocean circulation and ecosystems within the Chukchi Sea.

The spring phytoplankton bloom plays a major role in pelagic ecosystems; however, its dynamics is overlooked due to insufficient, highly-resolved observational data. Here we investigate the start, peak and decline of a two-week phytoplankton spring bloom in Frohavet, located at the coast of mid-Norway. We used observations from an uncrewed surface vehicle (USV) combined with buoy measurements, satellite images, discrete water sampling and modelling approaches. The spring bloom (March-June 2022) consisted of multiple peaks (up to 5 mg m-3), with a long peak in April, coincident with the period when the USV captured the temporal and spatial dynamics of the bloom. Short-term (5 days) episode of calm weather in the spring, such as clear skies and consistent low wind speed (< 7 m s-1) shoaled the mixed layer depth (< 15 m), after strong wind speed (average wind speed up to 20 m s-1 in March) and mixing events in winter. These rapid changes in the environment promoted the rapid development of the spring bloom - from 1 to 5 mg m-3 in 5 days. Likewise, the collapse of the bloom was rather quick, 1-2 days and coincides with low nitrate values and rapid increase in wind speed (> 10 m s-1), suggesting strong influence of the environment on phytoplankton dynamics during early stages of the spring bloom. Understanding the dynamics of the spring bloom is crucial for the management of marine resources. Integration of distinct observational platforms has the potential to unveil the environmental factors underlying phytoplankton bloom dynamics.

The Beaufort Sea has experienced a significant decline in sea ice, with thinner first-year ice replacing thicker multi-year ice. This transition makes the ice cover weaker and more mobile, making it more vulnerable to breakup during winter. Using a coupled ocean-sea-ice model, we investigated the impact of these changes on sea-ice breakup events and lead formation from 2000 to 2018. The simulation shows an increasing trend in the Beaufort Sea lead area fraction during winter, with a pronounced transition around 2007. A high lead area fraction in winter promotes a significant growth of new, thin ice within the Beaufort region while also leading to enhanced sea ice transport out of the area. The export offsets ice growth, resulting in negative volume anomalies and preconditioning a thinner and weaker ice pack at the end of the cool season. Our results indicate that large breakup events may become more frequent as the sea-ice cover thins and that such events only became common after 2007. This result highlights the need to represent these processes in global-scale climate models to improve projections of the Arctic.

Observations indicate that symmetric instability is active in the East Greenland Current during strong northerly wind events. Theoretical considerations suggest that baroclinic instability may also be enhanced during these events. An ensemble of idealised numerical ocean models, forced with northerly winds show that the short time-scale response (from two to four weeks) to the increased baroclinicity of the flow is the excitation of symmetric instability, which sets the potential vorticity of the flow to zero. The high latitude of the current means that the zero potential vorticity state has low stratification, and symmetric instability destratifies the water column. On longer time scales (greater than four weeks), baroclinic instability is excited and the associated slumping of isopycnals restratifies the water column. Eddy-resolving models that fail to resolve the submesoscale should consider using submesoscale parameterisations to prevent the formation of overly stratified frontal systems following down-front wind events. The mixed layer in the current deepens at a rate proportional to the square root of the time-integrated wind stress. Peak water mass transformation rates vary linearly with the time-integrated wind stress. The duration of a wind event leads to a saturation of mixing rates which means increasing the peak wind stress in an event leads to no extra mixing. Using ERA5 reanalysis data we estimate that between 1.5Sv and 1.8Sv of East Greenland Coastal Current Waters are produced by mixing with lighter surface waters during wintertime by down-front wind events. Similar amounts of East Greenland-Irminger Current water are produced at a slower rate.

Ensemble Kalman Filters (EnKFs), which assimilate observations based on statistics derived from samples of ocean states called ensemble, have become the norm for ocean data assimilation (DA) and forecasting. These schemes are commonly implemented with inflation and localization techniques to increase their ensemble spread and to filter out spurious long-range correlations resulting from the limited-size ensembles imposed by computational burden constraints. Such ad hoc methods were found not necessary in ensemble DA experiments with simplified ocean/atmospheric models and large ensembles. Here, we conduct a series of 1-year-long ensemble experiments with a fully realistic EnKF-DA system in the Red Sea using tens-to-thousands of ensemble members. The system assimilates satellite and in-situ observations and accounts for model uncertainties by integrating a 4km-resolution ocean model with ECMWF atmospheric ensemble fields, perturbed internal physics and initial conditions for forecasting. Our results indicate that accounting for model uncertainties is more beneficial than simply increasing the ensemble size, with the improvements due to large ensemble leveling off at about 250 members. Besides, and in contrast to what is commonly observed with simplified models, the investigated ensemble DA system still required localization even when implemented with thousands of members. These findings are explained by (i) amplified spurious long-range correlations produced by the low-rank nature of the ECMWF atmospheric forcing ensemble, and (ii) non-Gaussianity generated by the perturbed internal physical parameterization schemes. Large ensemble forcing fields and non-Gaussian DA methods might be needed to take full benefits from large ensembles in ocean DA.

Given the possibility of irreversible changes to the Earth system, technological interventions such as solar radiation management (SRM) are sometimes framed as possible climate emergency brakes. However, little knowledge exists on the efficacy of such disruptive interventions. To fill in this gap, we perform Community Earth System Model 2 (CESM 2) simulations of a SSP5-8.5 scenario on which we impose either gradual early-century SRM to stabilise surface temperatures or a rapid late-century cooling, both realised via stratospheric aerosol injection (SAI). While both scenarios cool Earth’s surface, we find that ocean conditions differ drastically. The rapid-cooling scenario fails to dissipate sub-surface ocean heat content (OHC), ends up in a weaker AMOC state and does not restore an ailing North Atlantic deep convection. Furthermore, the weakened AMOC state mediates the climate response to rapid SAI, thus inducing an interhemispheric temperature asymmetry. Our results advise caution when considering SAI as an emergency intervention.

Aotearoa New Zealand’s marine environment is heavily impacted by El Niño Southern Oscillation (ENSO). Little is known about the effect of ENSO on oceanographic properties in the area or their consequences on the distribution of marine organisms. Here we characterise the spatio-temporal variability of fine-scale fronts (< 10 km) in the area of Tīpaka Moana Te Moana Nui a Toi Hauraki Gulf (HG) and investigate the impact of dominant wind direction, seasonality, and ENSO phase. We processed Ocean-Land Color Instrument (OLCI) images from satellites between 2016-2022 with a fit-for-purpose version of the Belkin and O’Reilly Algorithm, specifically for frontal detection. We find coherent shifts in the position of fine-scale features, with ENSO phases alternatively separate and connect different parts of the Gulf during El Niño and La Niña respectively. Overall, fronts tend to co-locate with the 70 m and 40 m isobaths in the outer and inner HG respectively and their locations shift close or away from shore in response to changes in dominant wind direction. Furthermore, offshore fronts occurrences increase during winter and spring, and nearshore ones increase during summer and autumn. Our results sketch a first assessment of the distribution of fine-scale features that are likely to impact the distribution of important areas for pollutants dispersion and feeding areas for marine megafauna.

A document by Chao Gao. Click on the document to view its contents.

Ocean warming is a key factor impacting future changes in climate. Here we investigate vertical structure changes in globally averaged ocean heat content (OHC) in high- (HR) and low-resolution (LR) future climate simulations with the Community Earth System Model (CESM). Compared with observation-based estimates, the simulated OHC anomalies in the upper 700 m and 2000 m during 1960-2020 are more realistic in CESM-HR than -LR. Under RCP8.5 scenario, the net surface heat into the ocean is very similar in CESM-HR and -LR However, CESM-HR has a larger increase in OHC in the upper 250 m compared to CESM-LR, but a smaller increase below 250 m. This difference can be traced to differences in eddy-induced vertical heat transport between CESM-HR and -LR in the historical period. Moreover, our results suggest that with the same heat input, upper-ocean warming is likely to be underestimated by most non-eddy-resolving climate models.