Introduction

The Earth System is inherently coupled, but ocean heat uptake determines the Earth’s energy budget [1] and global sea level rise [2]. The ocean has a key role linking to other parts of the Earth System, e.g. sea surface temperatures (SSTs) affect atmospheric circulation and precipitation [3], ocean circulation [4] and biogeochemistry [5] determine the flux of carbon between the oceanic and atmospheric reservoirs, and upper ocean temperatures and circulation influence sea ice properties and dynamics, while sea ice also feeds back onto the ocean circulation via insulation and freshwater input. Global and regional sea levels are strongly influenced by high-latitude ocean processes associated with the basal melting of Antarctic ice shelves [6] and ocean-driven melting of the Greenland ice sheet [7].

Since the first review in 2000 [8], the ocean components of Earth System models (ESMs) have evolved considerably following the different phases of coupled model intercomparison projects, CMIP3 to CMIP5, and most recently CMIP6 [9]. Conservation of heat and salt, exact computation of sea level, and the improvement of water mass properties have been the main objectives that have led modelling groups to abandon rigid-lid formulations and time-invariant vertical grids, and to develop higher-order advection schemes. Community modelling platforms such as NEMO [10], MOM (https://mom-ocean.github.io/), or POP [11] have been extended to offer a larger choice of grids, numerical schemes, and parameterisations. Two new ocean models with flexible grids have been introduced in CMIP6: FESOM [12] and MPAS [13], and although CMIP6 models generally use quasi-geopotential grids (‘z-coordinates’), one model has terrain-following coordinates (INM-CM5) [14] and two models use Lagrangian hybrid coordinates following isopycnals in the interior: MOM6 [15] and Nor-ESM [16].

This paper reviews the state of knowledge of the importance of ocean resolution drawing on recent results from CMIP6 models (including HighResMIP [17]) and building on previous reviews [18,19,20] to address to what extent ocean resolution introduces uncertainty into climate variability and projections of future climate. Although the emphasis of this review is on multi-decadal climate timescales, many of our conclusions also have relevance for initialised predictions with coupled forecast models [18, 21]. The ‘Resolution in Ocean Components of CMIP6 Earth System Models’ section reviews the current status of resolution in ocean components of CMIP6 Earth System models. The ‘Impact of Ocean Resolution on Mean State, Variability, and Future Projections’ section reviews the impact of ocean resolution on the mean state, variability, and future projections of key metrics. The ‘Links to Other Aspects of the Earth System’ section reviews linkages between ocean resolution and other components of the Earth System reviewed in this issue. The ‘Advances in Parameterising the Mesoscale for Future Earth System Models’ section discusses advances in parameterisation of the mesoscale which are applicable for future developments of ESMs. The ‘Summary’ section summarises the review.

Resolution in Ocean Components of CMIP6 Earth System Models

Ocean models for climate, as atmospheric models, evolve constantly towards higher resolution. This is driven, in part, by the need to better represent strong western boundary currents such as the Agulhas Current, the Gulf Stream, or the Kuroshio, which play a key role in transporting heat from the equator to the poles. These currents have a typical width of 100 km, which means that grid meshes of 10–20 km are necessary to represent their dynamics. Improved horizontal resolution also allows for a better simulation of key straits such as the Gibraltar or Denmark Strait and their role in the inter-basin exchanges. Estimating grid resolution as the square root of the surface area of the Earth divided by the total number of grid points, the average resolution was 133 km for CMIP3, 87 km for CMIP5, and 58 km for CMIP6 (giving a timescale for doubling of ocean resolution in models of approximately 10 years [22]). The decrease is driven both by the refinement of the ESM grids and by the participation of more ocean general circulation models with very fine grids, in the 10–40-km range. This move has been strongly encouraged by HighResMIP [17], a protocol especially defined for high-resolution models. The OMIP, Ocean Model Intercomparison Project [23, 24], is also designed to help evaluate the impact of resolution on ocean simulations [25].

ESM development has to prioritise between allocating computational resources either to enhance resolution or to increase complexity (as well as considering the length of the spin-up and ensemble size). Although ESMs in CMIP6 have evolved in terms of their complexity [26,27,28,29], it is clear that most have retained ocean resolutions which require parameterisation of mesoscale eddies [30] as they fail to resolve the Rossby radius [17, 30]. The Rossby radius length scale is key to the representation of mesoscale eddies, boundary currents and fronts, and topographic flows [15, 18]. The length scales of the fastest-growing modes in the ocean have been shown to vary less strongly with latitude than suggested by the Rossby radius [31] with larger length scales observed in eastward-flowing currents. Nevertheless, taking the Rossby radius as a guide, three regimes of models [32] can be defined: eddy-parameterising (ocean resolution of O (50–100 km) and the current status for the majority of ESMs which mostly employ the Gent–McWilliams (GM) [33] parameterisation of mesoscale eddies), eddy-present (ocean resolution of O (25 km) with HighResMIP experiments at this resolution and some models of this resolution in the CMIP6 DECK experiment), and eddy-rich (ocean resolution of O (10 km) with some HighResMIP experiments at this resolution but noting that even models at this resolution do not resolve mesoscale eddies poleward of approximately 50°). Given that a major constraint for running at higher resolutions is the computational cost, some models also exploit variable resolution or nested high-resolution grids [34, 35].

A particular issue for models in the eddy-present regime is that the Rossby radius is resolved at lower latitudes but not at higher latitudes (or in shallow/shelf regions). This resolution is often referred to as the ‘grey zone’ when the Rossby radius is only marginally resolved and there is a question as to whether a parameterisation of mesoscale eddies should be included at this resolution. In CMIP6 models, the approach at the eddy-present resolution varies between models. For example, at 1/4° resolution, different model families make different choices about the use of GM [36, 37]. This is also an issue in eddy-rich models but at higher latitudes than in eddy-present models.

At even smaller spatial (and higher temporal) scales (100 m–10 km, hours to days), the parameterisation of submesoscale phenomena is also starting to be addressed in ESMs; the Fox-Kemper [38, 39] parameterisation of ocean surface boundary layer restratification by mixed layer instabilities has been followed by variants with different assumptions or for alternative submesoscale instabilities [40, 41] and implementations in ESMs [38, 42]. Submesoscale eddies need much higher resolution to resolve the smaller deformation radius of the mixed layer [43], but as demonstrated for the Agulhas system, there is evidence that the explicit simulation of submesoscale processes enhances the mesoscales [44]. Resolving submesoscale eddies is a target for climate simulations in future decades.

Impact of Ocean Resolution on Mean State, Variability and Future Projections

Western Boundary Currents

In typical CMIP5/6 models, with resolutions in the eddy-parameterising regime, simulated western boundary currents (WBCs) are much weaker and wider than observed, and the most significant biases in surface temperature are usually associated with incorrect presentation of WBCs in coupled models. Eddy-present models lead to a significant improvement in the representation of the WBCs as the ocean model becomes less diffusive/viscous with enhanced resolution, generally improving the simulation of the strength and position of WBCs such as the Gulf Stream, Kuroshio Current, and Agulhas Current [17, 25, 45,46,47,48,49,50]. For example, eddy-present and eddy-rich models show much better representation of the North Pacific subtropical gyre currents than eddy-parameterising models [18, 19, 25, 32, 49] where the Kuroshio separation and extension are often located at least 300–400 km further northward than observations. Despite general improvements in WBC representation due to model resolution, there remains dependence on the models’ numerics [51].

Mesoscale eddy activity is observed to be high in WBC regions, and the magnitude of eddy kinetic energy in these regions is strongly related to model resolution. Eddy tracking [52] demonstrates that an increasing number of eddies are simulated in WBC regions as resolution increases to the eddy-present and eddy-rich regimes. The magnitude of eddy kinetic energy in WBC regions is found to be strongly related to model resolution [37, 50, 51]. As mesoscale eddies potentially play a role in determining the strength of the gyre circulations and their low-frequency variability [53, 54], it is expected that decadal variability and sensitivity of the circulation to changes in wind stress will differ between eddy-rich and eddy-parameterised models.

In the case of the Agulhas Current where the leakage directly influences the hydrography in the Atlantic [55, 56], eddy-present and eddy-parameterising models overestimate the leakage by a factor of 2–3 [57]. There is also an indication that eddy-rich resolution is important for the determination of the relative role of warm and saline Agulhas leakage versus cold and fresh waters entering the Atlantic through the Drake Passage [58].

Some stand-alone high-resolution oceanic models tend to overestimate mesoscale eddies and underestimate the strength of WBCs, which is partly attributed to a lack of mesoscale air–sea interaction [59, 60] and partly attributed to incorrect representation of submesoscale motions in these models [61]. Eddy-present/rich coupled ocean–atmosphere models have substantially weaker eddies than eddy-present/rich ocean-only models without the surface current–wind stress interaction process because both the wind work and eddy-induced linear Ekman pumping dampen eddy kinetic energy when the surface currents are involved in the bulk formula [60, 62]. The Agulhas eddies’ pathway is also quite sensitive to the wind stress (‘relative’ versus ‘absolute’, depending on whether the ocean velocity is considered for the wind stress calculation) in forced simulations [51].

Beyond the eddy-rich resolution, further improvements (i.e. penetration) are found when the resolution is refined to ~ 1 km [50, 63, 64] and when a net release of available potential energy associated with mixed layer instability is responsible for the emergence of submesoscale eddies in wintertime [63]. At that resolution, submesoscale turbulence affects the kinetic energy exchanges and provides an inverse cascade of energy to larger scales [44, 64]. These extremely high–resolution experiments can be used in conjunction with machine learning techniques (e.g. deep learning) to design ocean eddy parameterisation for implementation in the coarser ocean models [65, 66] (see the ‘Advances in Parameterising the Mesoscale for Future Earth System Models’ section for further discussion).

Atlantic Meridional Overturning Circulation

The Atlantic meridional overturning circulation (AMOC) transports warm, buoyant water polewards in the Atlantic. Cooling at high latitudes, coupled with the relatively high salinity of the Atlantic, means that fluid parcels become sufficiently dense to convect and return equatorward at depth. This circulation drives a northward heat transport of about 1.2 PW at 26.5° N. Observations from the RAPID-MOCHA array since 2004 [67] show variability on all timescales with interannual variability typically of the order several Sv, with a larger temporary decrease of 4.7 Sv over 2009/2010 [68]. Most CMIP5 models [69] underestimate interannual variability and rarely if ever simulate such a large annual drop, though their daily variability can agree with observations [70, 71]. More recent eddy-present forced-ocean simulations [72] better capture the magnitude of interannual variability, suggesting that either ocean model resolution or atmospheric forcing/resolution is key. In CMIP6 HighResMIP and OMIP simulations (see Fig. 1a, b), the AMOC transport more often than not becomes stronger at higher ocean resolution [20], and in coupled models, this tends to be driven by enhanced convection in the Labrador Sea [73, 75]. For both CMIP5 models [76] and CMIP6 HighResMIP [73], there tend to be too deep mixed layers in the Labrador Sea in order for the AMOC strength to be comparable with observations (noting that this is particularly true in NEMO models, suggesting that the model structure as well as resolution may be a factor). When projecting future climate, this convection tends to reduce more quickly than that in the Nordic Seas [73], which means that the higher-resolution models have a stronger AMOC decline in the future. The effect of convection changes in the Nordic Seas is linked to the AMOC via the overflows.

Fig. 1
figure 1

Comparison of multi-model large-scale ocean metrics with observations. a AMOC and b northward heat transport (NHT) in the Atlantic at 26.5° N. c ACC transport (Sv). d An illustration of the impact of ocean model resolution on the ACC-AMOC simulation. Model simulations are from CMIP6 OMIP experiments [25] and from CMIP6 HighResMIP control-1950 experiments [73] (plus additional data not yet archived, see ‘Model Data’ for details). AMOC/NHT observations [67, 68] are shaded with the annual mean range between 2004 and 2017. ACC observations from the cDrake array [74] are shaded indicating 173 ± 11 Sv

The overflows of dense waters over the Greenland–Scotland Ridge (GSR), through Denmark Strait and the Iceland–Faroe channel, drive two-thirds of the AMOC [77]. Their dynamics are governed by small-scale physical processes such as ageostrophic flows in the bottom boundary layer, instabilities, and entrainment that are difficult to represent in numerical models [78, 79], especially at low resolution [80]. The relationship between the strength of the overflows and the AMOC is complex in models, because the AMOC depends on a number of processes [81]. Rather than the strength of the AMOC, its vertical structure critically depends on the depth to which overflows sink [80]. Also, the overflows tend to be stable over decadal timescales, which results in the AMOC decadal and subdecadal variabilities being mostly dependent on the water mass formation in the subpolar gyre [81]. On the other hand, the influence of the overflows is strong on the water mass properties and on the deep stratification of the subpolar gyre, which in turn influences the depth of wintertime convection [82].

The choice of vertical grid is expected to influence the representation of dense overflows. CMIP6 models generally employ z-coordinates (with variants such as time dependence of thickness to follow the motions of the free surface (z*) and partial cell topography in the deepest layer), but terrain-following coordinates and Lagrangian hybrid coordinates are also used [14,15,16]. Hybrid coordinates are designed to avoid spurious diapycnal mixing and better represent the entrainment downstream of the overflows. However, even in these models, the degree of improvement for the depth of the overflows and the AMOC vertical structure depends on the horizontal resolution and other model choices [15, 16, 81]. In all vertical coordinate systems, moving to higher horizontal resolution improves the representation of the overflows [16, 83] and their influence on AMOC, both because the bathymetry is more realistic and because the bathymetry is more realistic and because the dynamics of the deep boundary currents is better represented [84]. Higher vertical resolution, however, can degrade the solution in z-coordinate models, because it can increase spurious diapycnal mixing in downslope flows [83]. A complete analysis of overflows in HighResMIP ocean models is not yet available for a full assessment of the progress made relative to CMIP5.

Antarctic Circumpolar Current and Associated Southern Ocean Dynamics

Different aspects of the regional dynamics in the Southern Ocean are closely dynamically linked including the Antarctic circumpolar current (ACC), the Southern hemisphere upper and lower MOC cells, temperature and salinity structure and associated meridional gradients, and sea ice processes [85]. For example, the Southern Ocean mean state and its response to expected future changes in wind stress forcing, including the upper cell of its overturning circulation which influences heat and carbon uptake and isopycnal slopes that drive the ACC, result from a subtle balance between larger opposing wind-forced and eddy-driven cells [86]. Perhaps in consequence, representation of all of these ocean processes shows considerable sensitivities to ocean model resolution and configuration, even including fairly subtle changes to closure schemes which represent the impacts of sub–grid-scale processes [87, 88].

The simulation of the ACC is sensitive to model resolution in both forced and coupled simulations (Fig. 1c). Relative to observations from the cDrake array [74], models tend to fall on the lower end of observational uncertainty or underestimate the net ACC volume transport [85]. Models in the eddy-parameterised or eddy-rich regime perform better than those in the eddy-present regime which are particularly sensitive to the choice of uncertain grid-scale closures as seen in forced-ocean experiments [15, 36, 37].

The few CMIP6 models in the eddy-present regime provide a much more realistic representation of the frontal jets of the ACC [85] in contrast to the majority of models which have lower resolution and completely fail to represent any distinct fronts. These higher-resolution models, however, also show distinct counter flows (presumably linked to stationary eddies in Drake Passage) in multi-decadal means which are not evident in the observational means over similar periods (although they may be evident in sections from individual cruises). Increased model resolution allows for greater topographic detail to be resolved leading to improved topographic control of the ACC. In CMIP5, a weaker relationship between the ACC strength and position with westerly winds compared with CMIP3 was attributed to this control [89]. Further work is required to investigate this in CMIP6 and eddy-rich models.

There is much ongoing effort to better understand the strong dependence of the ACC on ocean resolution. Preliminary unpublished findings suggest that in Hadley Centre models, ACC transports using ocean resolution of 1/4° are rather sensitive to subtle changes in the configuration of grid-scale closures, such as viscosity, the use of weak GM, the use of localised partial slip at topographic boundaries, and even changes to bed stress. Analyses of perturbed parameter ensemble with a 1/4° ocean model [90], in which only atmospheric parameters are perturbed, also suggest that there may be a link between large-scale Southern Ocean SST biases; sea ice concentrations and associated meridional freshwater transports; near coastal salinity biases that drive a strong westward counter flow around the continental margin through the Drake Passage, weakening the total ACC transport; and Weddell Sea sea ice and salinity biases that cause a spurious polynya to form causing deep convection and large changes to Antarctic bottom water properties.

Several other studies have also noted multi-decadal temporal variations in ACC transports, of 10 to 30 Sv, linked to the occurrence of unrealistic large open-ocean polynya events, in both the Ross and Weddell Sea, which alter the density structure of the Southern Ocean through intense spurious open-ocean convection [85]. In 1/4° and 1/2° GFDL models [91, 92], super-polynyas in the Ross Sea have been shown to drive large centennial-scale variability in the Southern Hemisphere climate. These events are found in both the eddy-parameterising models used in CMIP6 and higher-resolution models [91,92,93]. Whether there is a resolution dependence on their frequency and characteristics is a topic of future study. With the exception of the 1974–1976 [94, 95] and 2016–2017 [96] polynya events observed in the Weddell Sea, there is no observational evidence to support the frequent open-ocean polynyas and associated intense open-ocean convection that are found in model simulations. Furthermore, there is no observational evidence of large open-ocean polynyas in the Ross Sea, where these events are commonly simulated in models. For these reasons, the dynamics associated with these events is a topic of focus as a means for improving simulations in future model development.

Substantial changes to the strength and latitudinal location of the westerly winds over the Southern Ocean have occurred and are expected to continue through a combination of ozone and greenhouse gas climate forcing. Future wind-driven changes in the strength of the ACC under climate change, and associated changes to the upper cell of the Southern Ocean MOC, may also depend strongly on ocean resolution. Idealised equilibrated ocean-only experiments suggest that at equilibrium in eddy-rich models, changes to the Ekman-driven overturning cell due to increased wind stress forcing are opposed by subsequent changes to the counter-rotating eddy-driven cell caused by increases in eddy kinetic energy [97] although the transient adjustment could be different to the equilibrium. This considerably reduces the sensitivity to changes in wind stress magnitude of both the upper cell of the Southern Ocean overturning circulation and, through its impact on isopycnal slopes, the strength of the ACC, processes termed ‘eddy compensation’ and ‘eddy saturation’, respectively. Evidence suggests that although the eddy response is not correctly represented in eddy-parameterising models, partial eddy compensation and eddy saturation are obtained in some models in which GM is allowed to vary spatially and temporally [87] and even to some extent with constant GM [98].

One interesting finding across the CMIP6 and associated HighResMIP ensemble is that many models with a (good) high ACC transport often have a (poor) low AMOC strength, and vice versa (Fig. 1d, b). This requires further investigation, but such a relationship could arise through either opposing responses to common grid-scale numerics (for example, more damping could result in weaker Gulf Stream transports but stronger ACC transports) or opposing dynamical links (for example, the Gnanadesikan model [99] suggests that a stronger AMOC should act to reduce the density gradients across the ACC, weakening its transport).

Sea Surface Temperature

SST is a key metric since it is the primary way that the ocean impacts the atmosphere and vice versa. The large-scale pattern of SST in ESMs demonstrates errors that are broadly similar across a range of models with cooling in the North Pacific, warming in the Southern Ocean, warming in the Eastern Boundary Upwelling zones, and a patch of very cold water (‘blue spot’) in the North Atlantic [26, 28, 29, 91, 100]. In many models, the magnitude of the Southern Ocean warming has been improved significantly through focussed effort largely to improve cloud biases [101]. Enhancing ocean resolution often leads to a reduction in the North Atlantic cold bias [102, 103] associated with the improvement in the position of the North Atlantic current, but forced-ocean experiments suggest that this is not always the case [15, 25, 36, 37]. In coupled models with eddy-present or eddy-rich resolutions [35, 49, 104, 105], there are typical but not uniform improvements in the Atlantic in both the tropics (both cold bias and zonal gradient) and mid-latitudes, the cold tongue in the tropical Pacific and the warm biases in the upwelling/stratocumulus regions off the western coasts of South America and Southern Africa, while the Southern Ocean surface warm bias tends to be increased.

Key impacts of parameterisations of submesoscale restratification in ESMs are shallower mixed layers particularly in wintertime, affecting air–sea fluxes, energy transfers, sea ice, and biogeochemistry [38, 106,107,108]. Without further alteration to other aspects of the ESM, most implementations see improved extratropical winter hemisphere mixed layers but degraded tropical and austral summer Southern Ocean biases [38, 42]. However, the use of the submesoscale parameterisation is not a good predictor as to whether models will have a good representation of the mixed layer [75, 109], which is perhaps unsurprising given the level of disagreement between models of boundary layer mixing [110].

Heat Uptake

Ocean and coupled models initialised from rest with climatological temperature and salinity generally gain or lose heat as they approach a quasi-equilibrium. The underlying drifts in temperature can be affected by both horizontal resolution [32, 104, 111] and the choice of vertical coordinate [15]. Models with horizontal resolutions in the eddy-present regime often have deep biases [32, 104, 111] that are worse than biases in models in either the eddy-parameterising or eddy-rich regime which may be linked to surface biases [101] propagated to the subsurface [111]. Excess interior mixing due to numerics [112] can be alleviated by moving from a depth coordinate to an isopycnal coordinate in the interior of the ocean [15], reducing heat uptake in the GFDL model by a factor of approximately 500 from 1 W/m2 when a hybrid depth–isopycnal vertical coordinate was used.

Work with simplified models [113, 114] suggests that heat uptake due to anthropogenic forcing will occur in the ventilated thermocline on the timescale of a few decades and in the deeper ocean on a longer timescale associated with changes in the AMOC and that insufficiently resolved eddies in the eddy-present regime could lead to reduced ocean heat uptake at depth. In idealised experiments, the Southern Ocean plays a particularly important role in global ocean heat (and anthropogenic carbon) uptake [114] (see also the ‘Antarctic Circumpolar Current and Associated Southern Ocean Dynamics’ section), and more realistic sensitivity experiments have suggested that Southern Ocean heat uptake is sensitive to horizontal resolution when moving from the eddy-parameterising to the eddy-present regime [115], to the near surface vertical resolution [116], and, in eddy-parameterising models, to the magnitude of the thickness diffusion parameter [117].

There are currently insufficient results from eddy-rich simulations to assess the impact of this resolution on projections of future ocean heat uptake. Projections of ocean heat uptake from eddy-rich models are likely to be limited by the length of the spin-up, as the magnitude of the underlying drift is likely to affect stratification and rates of ocean heat uptake [118]. A long spin-up which would be required for smaller drifts [119] is currently too expensive for eddy-rich models, but in the next decade, as DECK simulations are extended to include eddy-rich models, this will be an area for future research.

The uptake of anthropogenic CO2 over the historical period and into the twenty-first century is heavily influenced by the background ocean transport, which includes large-scale advection and eddy transport [120]. Like heat, a large fraction of the anthropogenic CO2 taken up by the oceans occurs in the Southern Ocean, and we can expect that the storage might be sensitive to horizontal resolution.

Links to Other Aspects of the Earth System

Sea Ice

Sea ice is comprised of floes and resembles a generally moving jumble of irregular, often interlocked, pieces of ice that vary in size from a few metres up to tens of kilometres. Crucially, sea ice is not a turbulent fluid, and so we would not expect performance to respond to grid resolution changes in the same way as in the ocean or atmosphere [121, 122].

Sea ice is highly non-Gaussian in nature and exhibits considerable heterogeneity in both time and space. To account for this, sea ice evolution in contemporary climate models is targeted at scales of ~ 100 km over periods of days to months and expressed in terms of local balances of conserved quantities such as mass and heat, with unresolved, small-scale processes handled using numerous parameterisations. This modelling framework, using the viscous-plastic (VP) model [123], or a derivative of it, is based upon an isotropic, plastic continuum approach whose validity relies upon statistical averages taken over a large number of floes [124]. Therefore, simply increasing the resolution of the sea ice model component will likely have little impact on the evolution of the sea ice per se (although when resolution is increased, parameterisations may be able to be better optimised). In spite of this, and of the fact that the continuum model assumptions break down at eddy-present resolutions, small-scale features can be obtained from VP continuum–based models through virtue of the high-resolution atmosphere and ocean boundary forcing. Simulations at kilometric resolutions using isotropic, plastic rheologies have been shown to generate linear kinematic features, such as leads and ridges, which look realistic, although are not well resolved and may not be oriented correctly [125].

Furthermore, being at the interface between the two, sea ice responds strongly to the forcing provided by the atmosphere and ocean components within climate models [126]. This means that changes in model resolution that lead to improvements in the mean state (i.e. bias reduction) or variability of the atmosphere or ocean models can lead to considerable improvements for the sea ice. For example, improvements to Southern Ocean circulation and SST biases can have a considerable impact on Antarctic sea ice cover [104, 111], while improved transport of warm Atlantic waters into the Arctic through Fram Strait can have a considerable impact on Arctic sea ice thickness [127].

Although continuum sea ice models are not expected to resolve small-scale dynamical features, they represent the heterogeneity of sea ice using a subgrid ice thickness distribution (ITD). Climate models with the most complexity in their sea ice components generally use a prognostic ITD to evolve the ITD explicitly at each model grid cell. This can be considered akin to a resolution increase for sea ice models and can have a considerable impact on heat exchanges over, and through, sea ice (Komuro and Suzuki, 2013). Inclusion of a prognostic ITD has been shown to have a considerable impact on sea ice evolution and feedback within climate models [128]. For example, enhancement of the (positive) ice–albedo feedback, coupled with suppression of the (negative) thickness–growth and thickness–strength (i.e. ridging) feedbacks, can lead to enhanced ice loss [129].

Biogeochemistry

Models of marine biogeochemistry are typically used within ESMs to represent the ocean’s role in the global carbon cycle and increasingly used to explore how climate change may impact marine ecosystems [130, 131]. Eddy-parameterising models still permit simulation of large-scale features and trends in marine biogeochemistry [132]. However, the ocean features that drive both the climate dynamics and the dynamics of many fish species ideally require much greater model resolution and regional realism [130, 133]. The computational cost of increased spatial resolution in models is compounded by the growing complexity of marine biogeochemistry models through expensive tracer advection [134]. Consequently, there are numerous ongoing activities to address biogeochemical complexity [135, 136] as well as the reduction of simulation cost [134, 137, 138]. These have important implications, for instance lengthening spin-up to reduce model drift [119] or in efficiently tuning model biogeochemistry [139]. Overall, model resolution and process complexity necessarily trade-off against simulation duration or experimental range. Choice of model sophistication may be clear where an end application requires only large-scale accuracy over long duration, or fine-scale accuracy over shorter periods. Many activities span this range, requiring an awareness in end users of the limitations that particular simulations impose in terms of spatial, physical, and biogeochemical accuracy. For example, results with idealised models suggest that the ocean carbon budget may exhibit reduced sensitivity to variations in Southern Ocean wind magnitude with explicit rather than parameterised eddies [140].

Ice Sheet Modelling

The horizontal size of Antarctic ice shelves ranges from the tiny Ferrigno ice shelf (~ 100 km2) to the giant Ross ice shelf (~ 500,000 km2) [141]. Basal melt can also be distributed over a range of length scales [142], concentrated in the first 20 km near the grounding line but with kilometre-scale variation due to the presence of a network of basal channels. Vertical length scales of the ice shelf also need to be considered: with buoyant plumes of O(1)m, simulated melt is very sensitive to the vertical resolution at the ice shelf base and vertical discretisation in the ocean model as well as the computation of the thermal driving [143]. Although many ocean models include ice shelf/ocean interactions [144, 145], current ESMs with eddy-parameterising and eddy-present resolutions are unable to capture all the relevant processes due to the horizontal and vertical length scales. One approach is to parameterise the circulation inside the ice shelf cavity (ice pump, ice shelf melt) on scales that the models are able to resolve [144, 146].

In current ESMs, representation of fjords around Greenland glaciers is missing. Explicit representation of processes at tide water glacier fronts requires very high horizontal and vertical resolution (O(1)m) [147] due to the typical size of the buoyant plumes generated by the subglacial runoff. Therefore, parameterisation of fjord dynamics and representation of ocean/tide water glacier interactions are needed to link the glacier front to the modelled ocean. Work with detailed models of fjord dynamics [148,149,150] (at resolutions much finer than ESMs) and analytical models [151, 152] are working towards parameterisations suitable for ESMs.

The impact of cold, fresh water resulting from ice sheet melting on the large-scale ocean circulation is also resolution dependent. Convection in the subpolar North Atlantic behaves differently in volumes and timing depending on whether melt water from Greenland glaciers finds its way into the interior of the Labrador Sea or is flushed within the Labrador Current towards the south [153, 154]. Since this boundary–interior exchange takes place via mesoscale eddies [155], it requires grid scales of 1/20°. In consequence, to correctly explore the potential response of the AMOC on ice sheet decay, it requires an ESM with eddy-rich resolution.

Advances in Parameterising the Mesoscale for Future Earth System Models

Current mesoscale eddy parameterisations in eddy-parameterising models, based on GM, mimic baroclinic instability and act as a net sink of available potential energy. As discussed in the ‘Impact of Ocean Resolution on Mean State, Variability, and Future Projections’ section, GM does not fully capture the physics of eddy saturation and compensation. It has been proposed that solving an explicit eddy energy budget is critical to understanding and correctly modelling eddy saturation [156]. The new GEOMETRIC eddy parameterisation follows such an approach, using the parameterised eddy energy to rescale GM [157] and looks promising in terms of reproducing both eddy saturation and eddy compensation [158].

As discussed in the ‘Resolution in Ocean Components of CMIP6 Earth System Models’ section, many models in the eddy-present regime do not incorporate GM but also do not explicitly resolve the mesoscale field, which can lead to less realistic behaviour in eddy-present models than eddy-parameterising ones. Scale-aware implementations of GM will allow the scheme to be used at eddy-present and eddy-rich resolutions without killing the eddy field [41]. Variability in eddy-present and eddy-rich models has also been shown to differ from that in eddy-parameterised models. This is to be expected as most parameterisations are designed to only capture the mean effects of eddies and not the variability of eddies, although some recent schemes attempt to also parameterise the variability [159,160,161,162]. The choice of resolving or parameterising the mesoscale at eddy-present (‘grey zone’) resolutions is likely to remain difficult.

There is currently no parameterisation implemented in climate models that mimics the important transfer of kinetic energy from small to large scales, namely kinetic energy backscatter, which affects the large-scale flow [163]. Recently, several studies using idealised numerical setups have developed energy backscatter parameterisations for ocean models. There are two main categories of parameterisations currently being developed:

  • A new set of momentum closures that have shown promise in mimicking kinetic energy backscatter. For example, stochastic eddy parametrisations have been developed for both eddy-parameterising and eddy-present simulations [160, 164]. The statistics of the stochastic models are crucial to mimicking eddy-mean flow interactions and improve the large-scale biases in ocean currents. Another example is a class of flow- and scale-aware parameterisations based on a non-Newtonian stress relation [161, 165, 166]. Finally, anti-viscous parameterisations have also shown improvements at eddy-permitting resolution [167, 168], when energetically constrained.

  • At eddy-parameterising resolution, part of the available potential energy extracted by the GM parameterisation can be re-introduced at a large scale. For example, studies have reinjected the energy lost by GM into the momentum equation using a simple anti-viscous term [169,170,171]. We note that this approach (a) could be combined with the new momentum closures which are more physically appropriate than an anti-viscous term [163] and (b) may have implications for computational cost due to a reduction in timestep required to satisfy the CFL criteria [171].

Beyond the mesoscale, considering parameterisation of the submesoscale for coarser resolution models will be key to future model improvement. In particular, parameterisation of other submesoscale impacts that are distinct from mixed layer eddy restratification, such as damping through submesoscale air–sea fluxes [172], lateral dispersion of pollutants and tracers [173], and submesoscale vertical transport below the mixed layer and nutrient pumping [174], is being developed. Different classes of submesoscale parameterisations, those specifically designed for use as subgrid schemes to carry forward cascades of energy, enstrophy, and tracer variance in mesoscale-resolving models, are in development [41, 161, 175, 176] and being prototyped in realistic simulations where they are shown to have benefits in less damping leading to better realism of energy and other dissipation [41, 177].

To satisfy the varying effective resolution versus a spatiotemporally variable mesoscale and submesoscale eddy scale, scale and flow awareness [30, 43] is essential ingredients of this class of parameterisations. One of the big advantages of scale-aware parameterisations is that they support the use of a hierarchy of resolutions as it avoids the need to retune parameterisations for each resolution [39, 41, 161].

Summary

The choice of ocean resolution in full Earth System models will always be limited by computational resources. Although there has been notable progress in increasing ocean resolution since CMIP3, the average resolution of the ocean component is still above 50 km with an effective resolution on the order of 300 km [178] which is more than five times greater than the resolution required to resolve the Rossby radius at mid-latitudes. Pioneering ocean simulations at eddy-rich resolution offers the potential to better evaluate errors and parameterisations in lower-resolution full Earth System models with the aim of better constraining their simulation of climate.

This review has demonstrated that both the mean state of the climate and the variability are sensitive to the choice of horizontal ocean resolution. The mean state is not uniformly improved by increased resolution, and the sensitivity is generally different across key metrics such as the Atlantic meridional overturning circulation and the Antarctic circumpolar current. This demonstrates that the numerical and parameterisation choices within configurations of ocean models remain important when producing the best possible representation of the ocean. In addition to horizontal resolution, vertical coordinates and resolution are also a key factor in ocean model performance, both in terms of capturing the baroclinic modes [179] and in particular regions determined by bathymetry such as overflows and ice cavities as well as in the surface boundary layer.

Particularly relevant to the use of ocean models as a component of Earth System models is that the choice of ocean resolution has effects beyond the ocean physics itself. Difficult choices with respect to resolution will need to be made in the future to satisfy the requirements for simulating the ocean biogeochemistry and capturing the details of ocean–ice sheet interactions as well as maintaining a computational cost that allows long multi-centennial simulations both for spin-up and projections. This is likely to require implementing variable horizontal resolution so that resolution can be placed in areas of high ocean eddy activity and critical frontal or topographic features.

Model Data

Model simulation output used in Fig. 1 can be obtained via the Earth System Grid Federation (ESGF) nodes for CMIP6 HighResMIP: HadGEM3-GC3.1 [180,181,182], ECMWF-IFS [182, 184], CNRM-CM6-1 [185, 186], CMCC-CM2-(V)HR4 [187, 188], EC-Earth3P [189, 190], MPI-ESM1-2 [191, 192], CESM1-3 [193, 194] and DOE-E3SM [195]. For the CMIP6 OMIP2, an archive of the model outputs and the scripts used to process the data is available at https://doi.org/10.5281/zenodo.3685918.