Impact of horizontal resolution on global ocean–sea ice model simulations based on the experimental protocols of the Ocean Model Intercomparison Project phase 2 (OMIP-2)

This paper presents global comparisons of fundamental global climate variables from a suite of four pairs of matched low- and high-resolution ocean and sea ice simulations that are obtained following the OMIP-2 protocol (Griffies et al., 2016) and integrated for one cycle (1958–2018) of the JRA55-do atmospheric state and runoff dataset (Tsujino et al., 2018). Our goal is to assess the robustness of climate-relevant improvements in ocean simulations (mean and variability) associated with moving from coarse (∼ 1^∘) to eddy-resolving (∼ 0.1^∘) horizontal resolutions. The models are diverse in their numerics and parameterizations, but each low-resolution and high-resolution pair of models is matched so as to isolate, to the extent possible, the effects of horizontal resolution. A variety of observational datasets are used to assess the fidelity of simulated temperature and salinity, sea surface height, kinetic energy, heat and volume transports, and sea ice distribution. This paper provides a crucial benchmark for future studies comparing and improving different schemes in any of the models used in this study or similar ones. The biases in the low-resolution simulations are familiar, and their gross features – position, strength, and variability of western boundary currents, equatorial currents, and the Antarctic Circumpolar Current – are significantly improved in the high-resolution models. However, despite the fact that the high-resolution models “resolve” most of these features, the improvements in temperature and salinity are inconsistent among the different model families, and some regions show increased bias over their low-resolution counterparts. Greatly enhanced horizontal resolution does not deliver unambiguous bias improvement in all regions for all models.

1 Introduction

A key decision in climate model design is the spatial resolution of different components. The global scope of the integrations and the centennial to millennial timescales and multiple scenarios required to capture changes in the climate set the problem; the spatial resolution is therefore the result of available computing. This decade, computing has become sufficiently powerful to make mesoscale-rich resolution affordable in ocean and sea ice models over most of the Earth, which allows the ocean model to simulate more intense internal variability than a lower-resolution model. As this new regime of coupled modeling is entered, it is important to understand both the behavior of ocean and sea ice models in a controlled framework and the benefits and challenges that come with higher resolution. This work introduces a set of matched numerical simulations in low- and high-resolution ocean and sea ice models with the latest forcing protocol. It is anticipated that these results will inform fully coupled modeling studies wherein ocean resolution varies and also that follow-on studies will build on our results to examine process and regional detail.

In 2016, an international group of ocean modelers behind the development and analysis of global ocean–sea ice models used as a component of the Earth system models in CMIP6 proposed an Ocean Model Intercomparison Project (OMIP) (Griffies et al., 2016). The essential element behind the OMIP is a common set of atmospheric and river runoff datasets for computing surface boundary fluxes to drive the ocean–sea ice models, many of which are used as components of coupled climate system models. The OMIP protocol is an outcome of the Coordinated Ocean–ice Reference Experiments (COREs), which assessed the performance of ocean–sea ice models (Griffies et al., 2009, 2014; Danabasoglu et al., 2014, 2016; Downes et al., 2015; Farneti et al., 2015; Wang et al., 2016a, b; Ilicak et al., 2016; Tseng et al., 2016; Rahaman et al., 2020) using the atmospheric and river runoff dataset of Large and Yeager (2009). However, this dataset has not been updated since 2009, and a new dataset (JRA55-do; Tsujino et al., 2018) has been developed for the OMIP based on the Japanese Reanalysis (JRA-55) product from Kobayashi et al. (2015) to ensure that it is regularly updated. This raw reanalysis product has been substantially adjusted to match reference states based on observations or the ensemble mean of other atmospheric reanalysis products as detailed in Tsujino et al. (2018) to create a suitable forcing dataset for ocean and sea ice models, referred to as JRA55-do. The continental river discharge is provided by a river-routing model forced by input runoff from the land surface component of JRA-55 adjusted to ensure similar long-term variabilities as in the CORE dataset (Suzuki et al., 2018). Runoff from ice sheets and glaciers from Greenland (Bamber et al., 2012, 2018) and Antarctica (Depoorter et al., 2013) are also incorporated. Tsujino et al. (2020) present an evaluation of the simulations from CMIP6-class global ocean–sea ice models forced with the JRA55-do datasets. This effort compares CORE-forced (i.e., OMIP-1) and JRA55-do-forced (i.e., OMIP-2) simulations considering metrics commonly used in the evaluation of global ocean–sea ice models to assess model biases.

Many features are very similar between OMIP-1 and OMIP-2 simulations, but Tsujino et al. (2020) identify many improvements in the simulated fields in transitioning from OMIP-1 to OMIP-2. For example, the sea surface temperature of the OMIP-2 simulations reproduces the global warming hiatus in the 2000s and the recent observed warming, which are absent in OMIP-1, partly because the OMIP-1 forcing stopped in 2009. The low bias in the sea ice area fraction in summer for both hemispheres in OMIP-1 is significantly improved in OMIP-2. The overall reproducibility of both seasonal and interannual variation in sea surface temperature and sea surface height (dynamic sea level) is also improved in OMIP-2. Tsujino et al. (2020) attribute many of the remaining model biases either to errors in representing important processes in ocean–sea ice models, some of which are expected to be mitigated by taking finer horizontal and/or vertical resolutions, or to shared biases in the atmospheric forcing. In this paper, we make a first attempt at quantifying the impacts of the models' horizontal resolution on biases.

Our goal is to assess the robustness of climate-relevant improvements in ocean simulations (mean and variability) associated with moving from coarse (∼ 1^∘) to eddy-resolving (∼ 0.1^∘) horizontal resolutions. Using the same atmospheric forcing (JRA55-do) for both low- and high-resolution configurations, we perform a multi-model analysis to identify the robust differences and improvements associated with increased resolution given the same forcing datasets. Within the ocean modeling community, it is usually assumed that high-resolution simulations should in general produce better results than low-resolution ones (Fox-Kemper et al., 2019). While this is clearly the case for surface currents and internal variability, we will show that greatly enhanced horizontal resolution does not necessarily deliver unambiguous bias improvement in temperature and salinity in all regions. It is important to note several caveats when interpreting the results presented in this paper: first, this is based on a limited number of numerical models (four), and second, because of the large computational cost associated with the high-resolution runs (factor of 1000 more expensive), only one JRA55-do cycle (1958–2018) is analyzed in this paper (versus six cycles for the coarse-resolution runs of Tsujino et al., 2020). Also, because of the short integration time, some of the results may not be robust (Atlantic meridional overturning circulation variability, deep-ocean circulation, etc.) (Danabasoglu et al., 2016). The layout of the paper is as follows. The models used in the comparison are described in Sect. 2. Section 3 highlights differences in the magnitude of the models' drift, while Sect. 4 focuses on the detrended interannual to decadal variability and the differences in the modeled ocean climates. The results are summarized and discussed in the final section.

Table 1 Model parameters for the low- and high-resolution configurations.

2 Description of the models

The CMIP6 OMIP-2 protocol does not include any specifications regarding model resolution, but most participating groups employ ocean models with horizontal resolution (∼ 1^∘) similar to what is used in the CMIP6 DECK experiments (Eyring et al., 2016) in order to achieve the required five-cycle spin-up (Tsujino et al., 2020). A high-resolution version of OMIP-2, with no multi-cycle spin-up requirement or well-defined protocols apart from the use of JRA55-do forcing, was informally organized by the CLIVAR Ocean Model Development Panel (OMDP) in 2019 to leverage the high-resolution (defined as being eddy-resolving over most of the globe, i.e., ∼ 1∕10^∘) work already being carried out by several of the modeling groups involved in the OMDP (coupled and uncoupled configurations). This study is an “intercomparison of opportunity” made possible by the handful of groups that were able to run parallel JRA55-do simulations at both eddy-parameterized (low) and eddy-resolving (high) resolutions. The high-resolution experiments are computationally expensive, and, when the call for comparison was made, each group leveraged known and proven configurations to perform the requested experiments. Furthermore, some groups had already completed the JRA55-do simulations at high resolution before this intercomparison was conceived. Given the large computational resources involved, rerunning those experiments to conform to a standard protocol was not an option. All experiments were configured using best practices, but each modeling group was empowered to choose what they thought was best and configured their high- and low-resolution configurations with similar parameters (see Table 1 for a detailed description of the parameters used in the low- and high-resolution model configurations). Ideally, only the horizontal resolution and associated physics should be changed to isolate the effects of horizontal resolution (Stewart et al., 2017), but this was not achievable for the present study since many of the low- and high-resolution simulations were often configured independently for distinct scientific goals and followed different development trajectories (e.g., vertical grids). It is also important to note that not all the models used the same climatology for the initial conditions, nor did they use the same wind stress formulation (absolute versus relative winds). When evaluating the drift of a numerical simulation, it is performed with respect to the climatology used to initialize the run.

2.1 FSU-HYCOM

FSU-HYCOM is a global configuration of the HYbrid Coordinate Ocean Model (HYCOM) (Bleck, 2002; Chassignet et al., 2003; Halliwell, 2004). The sea ice component is CICE version 4 (CICE4; Hunke and Lipscomb, 2010). The initial conditions in (potential) temperature and salinity are given by the Generalized Digital Environmental Model (GDEM4; Teague et al., 1990; Carnes et al., 2010). The Large and Yeager (2004) bulk formulation is used for turbulent air–sea fluxes except for the surface wind stress that is calculated without surface currents (absolute wind stress). No restoration is applied on the sea surface temperature. There is no parameterization of the overflows.

For the low-resolution configuration, FSU-HYCOM uses a tripolar Arakawa C grid of 0.72^∘ horizontal resolution with refinement to 0.33^∘ at the Equator (500 cells in the zonal direction and 382 in the meridional direction). The 2 min NAVO/Naval Research Laboratory DBDB2 dataset provides the bottom topography. A total of 41 hybrid coordinate layers are used, with potential density σ₂ target densities ranging from 17.00 to 37.42 kg m⁻³ (same configuration as Tsujino et al., 2020). The vertical discretization combines fixed pressure coordinates in the mixed layer and unstratified regions, isopycnic coordinates in the stratified open ocean, and terrain-following coordinates over shallow coastal regions. The surface salinity is restored to the monthly GDEM4 climatology over the entire domain with a salinity piston velocity of 50 m per 60 d, and the salinity flux at each time step is adjusted to ensure a net global flux of zero. Vertical mixing is the K-profile parameterization (KPP; Large et al., 1994) with a background diffusion of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹. Interface height smoothing by a biharmonic operator is used to correspond to Gent and McWilliams (1990), with a mixing coefficient determined by the grid spacing Δx (regular grid on a Mercator projection) times a velocity scale of 0.02 m s⁻¹ everywhere, except in the North Pacific and North Atlantic where a Laplacian operator with a velocity scale of 0.01 m s⁻¹ is used. Gent and McWilliams (1990) is not implemented where the FSU-HYCOM has coordinate surfaces aligned with constant pressure (mostly in the upper-ocean mixed layer) and there is no rotation of the lateral diffusion along neutral surfaces.

For the high-resolution configuration, FSU-HYCOM uses a tripolar Arakawa C grid of 0.08^∘ (1∕12^∘) horizontal resolution (4500 cells in the zonal direction and 3298 in the meridional direction). The model bathymetry is derived from the 30 arcsec GEBCO08 dataset. A vertical resolution of 36 hybrid layers, with σ₂ target densities ranging from 26.00 to 37.24 kg m⁻³, is used. The 36-layer high-resolution configuration was at the time our default configuration and was retained to compare to previous runs performed with other atmospheric forcing datasets. The 41-layer coarse-resolution runs were performed for inclusion in Tsujino et al. (2020) using the latest vertical grid with all the additional layers in the upper ocean. While the vertical resolution is lower in both configurations than recommended by Stewart et al. (2017) for z-coordinate models, the statistics of eddy scale and the vertical structure of the resolved eddy motions are well-captured with this layer discretization when compared to a z-coordinate model with 300 levels (Ajayi et al., 2020). The surface salinity is restored to the monthly GDEM4 climatology over the entire domain with a salinity piston velocity of 50 m per 60 d, and the salinity flux at each time step is adjusted to ensure a net global flux of zero. The KPP model (Large et al., 1994) with a background diffusivity of 10⁻⁵ m² s⁻¹ provides the vertical mixing. The horizontal advection uses a second-order flux-corrected transport scheme. An interface height smoothing is applied through a biharmonic operator (with a velocity scale of 0.015 m s⁻¹).

2.2 NCAR-POP

The NCAR-POP model is based on the ocean component of the Community Earth System Model version 2 (CESM2; Danabasoglu et al., 2020) and is a global configuration of the Parallel Ocean Program version 2 (POP2; Smith et al., 2010) with several modifications to the model physics and numerics, including improved treatment of continental freshwater discharge into unresolved estuaries (Sun et al., 2019) and a new parameterization of Langmuir mixing (Li et al., 2016). The sea ice component of CESM2 is CICE version 5.1.2 (CICE5; Hunke et al., 2015), which features new mushy-layer thermodynamics (Turner and Hunke, 2015) with prognostic sea ice salinity and an updated melt pond parameterization (Hunke et al., 2013). CICE5 uses the same horizontal mesh grid as the POP2 configuration to which it is coupled. The initial conditions are given by the World Ocean Atlas 2013 (WOA13; Locarnini et al., 2013; Zweng et al., 2013). The surface stress is a function of ocean surface velocity (relative wind stress), and sea surface salinity is restored to WOA13 monthly climatology with a piston velocity of 50 m over 1 year (after subtraction of the global mean restoring). Both configurations use a precipitation factor, computed once per year, to prevent salinity drift as discussed in Appendix C of Danabasoglu et al. (2014).

For the low-resolution configuration, NCAR-POP utilizes a dipole mesh grid with the grid northern pole displaced into Greenland. The horizontal resolution (nominal 1^∘) is uniform in the zonal direction (1.125^∘) but varies in the meridional direction (from 0.27^∘ at the Equator to ∼ 0.5^∘ at midlatitudes). The z-coordinate vertical grid has 60 levels, going from 10 m at the surface, to 250 m in the deep ocean, and to a maximum depth of 5500 m. The subgrid-scale closures and parameter settings used in this configuration are well-documented (Danabasoglu et al., 2012, 2014, 2020); some of the details are listed here to facilitate model intercomparison. The low-resolution NCAR-POP employs the skew-flux form of the GM isopycnal transport parameterization (Griffies, 1998), with depth-dependent thickness and isopycnal diffusivities (Ferreira et al., 2005; Danabasoglu and Marshall, 2007) from roughly 3000 m² s⁻¹ in the near surface to 300 m² s⁻¹ in the deep ocean. Near-surface mesoscale diabatic fluxes are also parameterized (Ferrari et al., 2008; Danabasoglu et al., 2008) with diffusivity set to 3000 m² s⁻¹, while the near-surface restratification effects of submesoscale mixed layer eddies are parameterized following Fox-Kemper et al. (2008, 2011). The K-profile vertical mixing scheme (KPP; Large et al., 1994) is used to prescribe the vertical viscosity and diffusivity coefficients as detailed in Danabasoglu et al. (2012). Their background values have a latitudinal structure but are constant in the vertical. Instead, a tidal mixing parameterization is used to include the effects of enhanced abyssal mixing from tidally generated breaking waves (Jayne, 2009). The low-resolution NCAR-POP (but not the high-resolution version) uses an overflow parameterization to represent the density-driven flows of the Denmark Strait, Faroe Bank Channel, the Ross Sea, and the Weddell Sea (Danabasoglu et al., 2010). This configuration uses hourly coupling rather than the daily coupling used in previous CORE publications (e.g., Danabasoglu et al., 2014).

For the high-resolution configuration, NCAR-POP utilizes a tripole grid with the grid northern poles in North America and Asia. It is based on versions documented in McClean et al. (2011) and Small et al. (2014) but has been updated to the CESM2 code base. The sea ice component, however, used CICE4 physics (i.e., excluding the CICE5 developments mentioned above) in order to maintain consistency with other high-resolution simulations that were run with earlier versions of CESM. The horizontal grid (nominal 0.1^∘) varies from 11 km at the Equator to 2.5 km at high latitudes, and the vertical grid (62 levels) is the same as that used in the low-resolution NCAR-POP but extends deeper to 6000 m. The additional two vertical levels (250 m thick) increase the maximum ocean depth from 5500 to 6000 m, allowing for a more realistic representation of deep-ocean features resolved by the 0.1^∘ grid. The high-resolution NCAR-POP uses a partial bottom cell formulation of the vertical grid for more accurate representation of bathymetry. In this configuration, biharmonic horizontal mixing is used for tracers and momentum, but there is otherwise no parameterization of eddy-induced mixing. New features in the CESM2 version include the use of half-hour coupling and a modified virtual salt flux formulation that uses a local reference salinity (the estuary parameterization used in the low-resolution version is not used in the high-resolution configuration, but the latter does use new methods for redistributing continental freshwater fluxes over several vertical layers near the surface). No Langmuir mixing parameterization is used in the high-resolution configuration.

2.3 AWI-FESOM

AWI-FESOM is a global configuration of the Finite Element/volumE Sea ice–Ocean Model (FESOM) version 1.4 (Wang et al., 2014; Danilov et al., 2015) of the Alfred Wegener Institute Climate Model (AWI-CM; Sidorenko et al., 2015, 2018; Rackow et al., 2018, 2019; Sein et al., 2018). Both the ocean and sea ice modules work on unstructured triangular meshes (Danilov et al., 2004; Wang et al., 2008), thus allowing for multi-resolution simulations. The tracer equations employ a flux-corrected advection scheme and the KPP scheme (Large et al., 1994) for vertical mixing. The background vertical diffusivity is latitude- and depth-dependent (Wang et al., 2014). Mesoscale eddies are parameterized by using along-isopycnal mixing (Redi, 1982) and Gent–McWilliams advection (Gent and McWilliams, 1990) with vertically varying diffusivity as implemented in Danabasoglu et al. (2008). The eddy parameterization is switched on when the first baroclinic Rossby radius is not resolved by local grid size. In the momentum equation, the Smagorinsky (1963) viscosity in a biharmonic form is applied. The sea ice thermodynamics follow Parkinson and Washington (1979), with a prognostic snow layer and snow to ice conversion. The Semtner (1976) zero-layer approach, assuming linear temperature profiles in both snow and sea ice, is used in this model version. The elastic–viscous–plastic (EVP; Hunke and Dukowicz, 1997) rheology is used with modifications that improve the convergence (Danilov et al., 2015; Wang et al., 2016c). The sea surface salinity (SSS) is restored to monthly PHC3 climatology with a piston velocity of 50 m over 300 d outside the Arctic Ocean and 3 times weaker inside the Arctic Ocean. The air–sea turbulence fluxes are calculated using the bulk formulation of Large and Yeager (2009). The full ocean surface velocity is used in the calculation of wind stress (relative wind stress). The initial conditions are derived from PHC3.0 (Steele et al., 2001).

Figure 1 Horizontal resolution (km) of the two FESOM grids: (a) low resolution and (b) high resolution.

For the low-resolution configuration, AWI-FESOM uses a nominal 1^∘ bulk horizontal resolution, which has been used in previous CORE-II simulations (e.g., Wang et al., 2016a, b), with the North Atlantic subpolar gyre region and Arctic Ocean set to 25 km (see Fig. 1a) and a nearly equatorial resolution of $\mathrm{1}/{\mathrm{3}}^{\circ }$.

For the high-resolution configuration, AWI-FESOM uses the grid introduced by Sein et al. (2016). As the variability of sea surface height (SSH) can manifest the variability of mesoscale eddies, the horizontal resolution is scaled with the strength of the observed SSH variability on this grid. In particular, the resolution is about 10 km along the western boundary currents, the Antarctic Circumpolar Current, and the Agulhas Current region (Fig. 1b). Before generating this grid, the field of SSH variance is smoothed spatially to make sure that the main currents are within high-resolution regions even if their positions change to some extent. The resolution is also increased along the coast and where the observed sea ice concentration variability is high. This multi-resolution grid has about 1.3 million surface grid points, similar to the size of a uniform $\mathrm{1}/{\mathrm{4}}^{\circ }$ mesh. In both setups, 46 z levels are used, with 10 m spacing in the upper 100 m. This is slightly less than the 50 levels recommended by Stewart et al. (2017) to resolve the first baroclinic mode in a z-coordinate model.

2.4 IAP-LICOM

IAP-LICOM is a global configuration of the LASG/IAP Climate system Ocean Model (LICOM) (Zhang and Liang, 1989; Liu et al., 2004, 2012; Yu et al., 2018; Lin et al., 2020) developed by the Institute of Atmospheric Physics (IAP) from the Laboratory of Atmospheric Sciences and Geophysical Fluid Dynamics (LASG) of the Chinese Academy of Sciences (CAS). LICOM is the ocean component of the Flexible Global Ocean–Atmosphere–Land System model (FGOALS; e.g., Li et al., 2013; Bao et al., 2013) and of the CAS Earth System Model (CAS-ESM; Minghua Zhang, personal communication, 2019). Version 3 of LICOM (LICOM3) is coupled to CICE4 through the NCAR flux coupler 7 (Craig et al., 2012; Lin et al., 2016). The surface salinity is restored to the monthly PHC3.0 climatology over the entire domain with a salinity piston velocity of 50 m per 4 years (50 m per 30 d for the sea ice regions). However, there is no freshwater flux normalization. The full ocean surface velocity is used in the calculation of wind stress (relative wind stress). The vertical viscosity and diffusion coefficients in the mixed layer are computed by the scheme of Canuto et al. (2001, 2002) with background values of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ and an upper limit of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{2}}$ m² s⁻¹. The tidal mixing scheme of St. Laurent et al. (2002) was recently implemented in LICOM3 by Yu et al. (2017). The initial conditions are derived from PHC3.0 for the low-resolution configuration and from the Mercator Ocean reanalysis for the high-resolution configuration.

For the low-resolution configuration, IAP-LICOM uses a Murray (1996) tripolar Arakawa B grid with two north “poles” at 65^∘ N, 65^∘ E and 65^∘ N, 115^∘ W with a resolution of approximately 1^∘ in both the zonal and meridional directions (360×218 grid points). The vertical grid uses the η coordinate (Mesinger and Janjic, 1985) with 30 levels. The horizontal viscosity consists of a Laplacian with a coefficient of 5400 m² s⁻¹. The mesoscale eddies are parameterized using the isopycnal tracer diffusion scheme of Redi (1982) and the eddy-induced tracer transport scheme of Gent and McWilliams (1990), with tapering factors as in Large et al. (1997) and the buoyancy-frequency-related (N²) thickness diffusivity of Ferreira et al. (2005).

Figure 2 Time evolution of the domain-averaged kinetic energy (10⁻⁴ m² s⁻²) for all experiments.

For the high-resolution configuration (Li et al., 2020), IAP-LICOM uses the same Murray (1996) tripolar grid as in the low-resolution version, but with two north “poles” at 55^∘ N, 95^∘ E and 55^∘ N, 85^∘ W and with a resolution of $\mathrm{1}/{\mathrm{10}}^{\circ }$ (11 km zonally and varying from 11 km at the Equator to 8 km in midlatitudes – 3600×2302 grid points). There are 55 levels in the vertical, which is just above the minimum recommended by Stewart et al. (2017) to resolve the first baroclinic mode. The additional 25 levels from the low-resolution configuration increase resolution in the deep ocean and improve the simulation of deep circulation and Atlantic meridional overturning circulation (AMOC) transport. The horizontal viscosity consists of a biharmonic operator with a coefficient of $-\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{10}}$ m⁴ s⁻¹. The Gent and McWilliams (1990) scheme is turned off and the tracers use a biharmonic isopycnal diffusivity with a coefficient of $-\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{10}}$ m⁴ s⁻¹. It important to note that only the thermodynamic part of CICE4 (no dynamics) was used in the high-resolution IAP-LICOM.

3 Temporal evolution and drift

3.1 Mean kinetic energy

Figure 2 shows the evolution of the domain-averaged mean kinetic energy for all experiments (solid lines for the high-resolution experiments, dashed lines for low-resolution experiments) from 1958 to 2018. The evolution is very similar between different models; all exhibit a quick spin-up of the kinetic energy in the first 5 years, which levels off for the rest of the integration. Not surprisingly, the total kinetic energy is significantly higher for the high-resolution experiments than the low-resolution experiments. For the high-resolution configurations, the FSU-HYCOM has the highest kinetic energy, with a globally averaged value of $\sim \mathrm{35}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻², and the IAP-LICOM has the lowest kinetic energy, with a globally averaged value of $\sim \mathrm{15}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻² in the high-resolution configuration. For comparison, a previous high-resolution $\mathrm{1}/{\mathrm{10}}^{\circ }$ global simulation, performed with an older version of POP by Maltrud and McClean (2005) using daily NCEP/NCAR reanalysis forcing and absolute winds in wind stress calculations, has a global averaged kinetic energy at 25–$\mathrm{30}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻² (see their Fig. 1). The higher kinetic energy in FSU-HYCOM can be partially explained by the wind stress formulation, which does not take into account the ocean current velocities (absolute winds), while the other three models do (relative winds). The latter has an “eddy-killing” effect that can reduce the total kinetic energy by as much as 30 % (see Renault et al., 2020, for a review). This is roughly the difference that is seen between FSU-HYCOM and NCAR-POP, and POP with absolute winds in the wind stress (Maltrud and McClean, 2005) has a level of kinetic energy that is close to FSU-HYCOM. But even with the highest resolution used here (∼ 0.1^∘), the total kinetic energy remains significantly lower that what can be inferred from observations and higher-resolution models (closer to $\mathrm{50}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻²; Chassignet and Xu, 2017). The increase in total kinetic energy from the low- to the high-resolution configuration is approximately a factor of 4 for all models, except for AWI-FESOM (factor of 2 only). This is probably because the high-resolution AWI-FESOM has a variable grid spacing (Fig. 1b) and does not resolve the Rossby radius of deformation everywhere. The level of total kinetic energy is substantially lower in IAP-LICOM, most likely because of the two-step shape preservation advection scheme used in the momentum equations and higher diffusivity (Table 1).

Figure 3 Time evolution from the initial conditions of the global steric sea surface height change for all experiments.

Figure 4 Time evolution of global temperature (^∘C) and salinity (psu) change (relative to initial conditions) for all experiments (depth-integrated, upper 700 m, upper 2000 m).

3.2 Global mean sea level, temperature, and salinity

As stated by Griffies et al. (2014), “[t]he CORE [and OMIP] protocols (Griffies et al., 2009, 2016; Danabasoglu et al., 2014) introduce a negligible change to the liquid ocean mass (non-Boussinesq) or volume (Boussinesq), and the salt should remain nearly constant (except for relatively small exchanges with the sea ice)”. Changes to the simulated global mean sea level should arise only because of thermosteric effects (i.e., changes in ocean heat content and redistribution of heat) in simulations that preserve salt content (i.e., that have zero net surface freshwater flux). In most models, the global sea level time evolution (Fig. 3) is dominated by changes in the global mean ocean temperature (Fig. 4a). IAP-LICOM is the exception, in which the global sea level shows a downward trend until ∼ 1990 and then slowly rises. This is due to an increase in global mean salinity (Fig. 4b), which dominates the global sea level changes despite a large increase in global volume mean temperature (Fig. 4a). This increase is explained by the lack of surface-restoring salt flux normalization in IAP-LICOM (see Sect. 2.4). For details on the salt flux normalization used in the other models, the reader is referred to Appendix B.3 of Griffies et al. (2009) and Appendix C in Danabasoglu et al. (2014), which describe the techniques used to ensure that there is no net salt added to or removed from the ocean–sea ice system (Griffies et al., 2014). One can also note that most models show an increase in global temperature and sea surface height after 1980–1990, associated with rising air temperature.

An increase in the horizontal resolution does not necessarily imply a reduction in temperature and salinity drift, and no coherent picture emerges from the comparison. If one focuses on the time evolution of the globally averaged temperature in the upper 700 m (Fig. 4c), there are only small changes for AWI-FESOM and IAP-LICOM, while FSU-HYCOM shows a significant reduction in the drift and NCAR-POP an increase as resolution is increased. For salinity in the upper 700 m (Fig. 4d), the increase in resolution significantly reduces the drift in NCAR-POP and FSU-HYCOM; there are no changes for IAP-LICOM and a significant increase for AWI-FESOM. It is important to note that, while the salt flux restoring may differ among the four models, it remains the same for each model as the resolution is increased. When considering the upper 2000 m (Fig. 4f), there is a significant reduction in the global salinity drift for NCAR-POP and FSU-HYCOM, less so for IAP-LICOM, and no changes for AWI-FESOM. We note that most high-resolution models (FSU-HYCOM being the exception) warm faster over the upper 2000 m and global temperature than their lower-resolution partners, which is not true for the upper 700 m. Figure 5 contains identical data as in Fig. 4, except that, by rebasing anomalies to the final year, it highlights the forced ocean variations of the last few decades of simulation by comparing the modeled global temperature and salinity as well as heat and salt content change to that of the World Ocean Atlas 18 (WOA18). The comparison shows that high resolution improves the fidelity and reduces the spread of forced ocean heat content change (particularly 0–2000 m heat content) over recent decades. Figure 6 shows in more detail the evolution of the global temperature and salinity from the initial conditions as a function of depth. AWI-FESOM shows the smallest changes in temperature throughout the water column but the largest in salinity despite having a salt flux adjustment to ensure that the global salinity remains constant (shown in Fig. 2). There is a significant freshening in the upper 400 m compensated for by a salinification in the deeper ocean. This is more pronounced in the high-resolution experiment. Increasing the resolution significantly improves the drift in both temperature and salinity in FSU-HYCOM but not so much in the other simulations. One could actually argue that the drift is stronger in NCAR-POP with a significant warming in the upper 400 m and freshening in the upper 100 m. While it is beyond the scope of this paper, additional insights could be gained by computing vertical heat and salt budgets as in Griffies et al. (2015) and von Storch et al. (2016). In the next section, we investigate in more detail the evolution of temperature and salinity as a function of depth and time by ocean basin.

Figure 5 Time evolution of global temperature (^∘C) and salinity (psu) change (relative to the year 2018) for all experiments and the World Ocean Atlas (WOA18) (depth-integrated, upper 700 m, upper 2000 m).

3.3 Temperature and salinity bias (depth vs. time) by ocean basin

Figures 7, 8, and 9 show the time evolutions of the horizontal mean depth profiles of temperature and salinity for the Atlantic, Pacific, and Indian Ocean, respectively. The color bar is the same in all figures, including Fig. 6 (global), therefore allowing for a qualitative estimate of where the drift is most significant. To a large extent, the time evolution in each of the major ocean basins mimics that of the global but with some significant differences. For the Atlantic Ocean (Fig. 7), there is a surface freshening in the upper 100 to 200 m as resolution is increased in FSU-HYCOM, NCAR-POP, and AWI-FESOM. IAP-LICOM, on the other hand, becomes more saline and warmer in the upper 100 m. The latter is accompanied by a large freshening and cooling below 100 m to approximately 600 m. Overall, the drift is smaller in the Pacific for most models, with a freshening in the upper ocean for all models as resolution is increased. The exception is again IAP-LICOM, which shows a significant increase in salinity in the upper 200 m, and this could be a consequence of the fact that there is no zero normalization of the surface-restoring salt flux. The temperature bias in the Pacific Ocean is similar to that of the global ocean. In the Indian Ocean (Fig. 9), FSU-HYCOM and NCAR-POP exhibit larger biases than AWI-FESOM and IAP-LICOM, but those biases are smaller at the higher resolution.

Figure 6 Time evolution of global temperature (^∘C) and salinity (psu) change (relative to initial conditions) as a function of depth for all experiments.

Figure 7 Time evolution of Atlantic Ocean temperature (^∘C) and salinity (psu) change (relative to initial conditions) as a function of depth for all experiments.

Figure 8 Time evolution of Pacific Ocean temperature (^∘C) and salinity (psu) change (relative to initial conditions) as a function of depth for all experiments.

Figure 9 Time evolution of Indian Ocean temperature (^∘C) and salinity (psu) change (relative to initial conditions) as a function of depth for all experiments.

3.4 Mean temperature and salinity bias

The drift plots from Sect. 3.2 and 3.3 indicate that temperature and salinity bias structures are well-established within the first 2 decades of the spin-up so that time averages computed over the latter decades of the simulations should provide a reasonable estimate of the stationary biases characterizing each model. Figure 10 shows latitude–longitude maps of surface temperature and salinity differences computed over the 1980–2018 time period with respect to WOA18. All models exhibit reductions in sea surface temperature (SST) bias as resolution is increased. The largest bias reductions are seen in the Southern Ocean, in the Agulhas retroflection region, and along the Gulf Stream extension in the North Atlantic. All these changes are presumably primarily related to an improved representation of the pathways of strong surface currents in those regions, which could result not only from differences between resolved and parameterized eddies but also from differences in resolved bathymetry. SST bias in eastern boundary upwelling regions has been shown to be strongly sensitive to atmospheric forcing resolution (Tsujino et al., 2020), somewhat sensitive to the regridding technique used for near-surface winds (Small et al., 2015), and also slightly sensitive to ocean resolution for seasonal means (Kurian et al., 2020). However, Fig. 10 shows a minimal change in annual mean SST bias in the upwelling zones off the west coasts of Namibia, Peru, and North America. In addition to the Gulf Stream, all models exhibit some degree of bias reduction in the Brazil–Malvinas Confluence zone at high resolution, but no such systematic improvement is evident in the Kuroshio–Oyashio extension region. Sea surface salinity bias is generally reduced globally in the high-resolution simulations, with the exception again being AWI-FESOM, which exhibits a more pronounced negative salinity bias in line with the enhanced near-surface salinity drift at high resolution in that model (see Figs. 6–9). The SSS bias in the Arctic results from the differences in salinity between WOA18 and the climatologies used in the salinity restoring (GDEM4, WOA13, or PHC3.0; see Sect. 2 for details).

Figure 10 Modeled surface temperature (^∘C) and salinity (psu) difference from WOA18 over the 1980–2018 time period.

The vertical structure of mean temperature and salinity bias is displayed as zonal averages across the different ocean basins in Figs. 11–13. This reveals that the near-surface global drift toward warming seen in FSU-HYCOM and NCAR-POP is partly related to a degradation of the tropical thermocline in all basins. One hypothesis is that the thermocline bias is related to the representation of vertical eddy heat flux (Griffies et al., 2015), which tends to be stronger and more realistic in high-resolution simulations (see Fig. 5). In FSU-HYCOM, higher resolution improves the situation, while in NCAR-POP, the tropical thermocline temperature bias actually gets worse with resolution. The degradation in thermocline bias in the POP high-resolution version could also be due to missing submesoscale physics, which are parameterized in the low-resolution configuration but absent in the high-resolution configuration. In both FSU-HYCOM and NCAR-POP, and to some extent IAP-LICOM, large extratropical and polar temperature biases associated with intermediate and deep waters decrease with enhanced resolution. In AWI-FESOM, both high- and low-latitude biases get worse. It turns out that the relatively low near-surface temperature drift seen in AWI-FESOM (Figs. 6–9) is related to large, but largely compensating, anomalies of opposite sign at each depth level. Cold–warm anomalies that develop in the tropics and at high latitudes in that model at both resolutions tend to disappear in global or basin-wide means. Such a compensation does not occur for salinity, which is too fresh in the thermocline and too saline in deeper intermediate waters in AWI-FESOM, a characteristic that gets worse as resolution is increased. IAP-LICOM exhibits some large changes in the sign and structure of zonal mean bias as resolution changes, but as with AWI-FESOM, it does not lend strong support to the hypothesis that model temperature and salinity bias can be reduced by increasing model horizontal resolution. As already mentioned, the NCAR-POP model exhibits improved representation of high-latitude intermediate and deep waters, but degraded representation of the tropical thermocline, as resolution is enhanced. FSU-HYCOM is the only model that shows a nearly ubiquitous reduction in temperature and salinity bias in all ocean basins in the high-resolution version; for the other models, greatly enhanced horizontal resolution does not deliver unambiguous bias improvement in all regions. This is a rather unexpected result.

Figure 11 Global zonal temperature (^∘C) and salinity (psu) difference with the climatology used to initialize the model.

Figure 12 Atlantic zonal temperature (^∘C) and salinity (psu) difference with the climatology used to initialize the model.

Figure 13 Indo-Pacific zonal temperature (^∘C) and salinity (psu) difference with the climatology used to initialize the model.

3.5 Northern and Southern Hemisphere sea ice volume and extent

Due to observational constraints, sea ice is often quantified and monitored by passive microwave satellites in terms of sea ice extent (SIE), defined as the sum of all grid areas with 15 % or higher sea ice concentration. Figure 14a–b display the modeled annual mean Northern and Southern Hemisphere sea ice extent for all simulations and include comparisons to the latest observations from the National Snow and Ice Data Center (NSIDC) (Fetterer et al., 2017). Similar to the observations of the 5th and 6th CMIPs (Stroeve et al., 2012; Notz et al., 2020; Shu et al., 2020), the multi-model mean is representative of the passive microwave mean SIE from 1978 to 2018, with large inter-model differences. All models show a clear decline in SIE in the Northern Hemisphere and weaker trends in the Southern Hemisphere. Because observed sea ice extent is highly correlated with near-surface air temperature (e.g., Olonscheck et al., 2019), the consistency of the temporal variability between different simulations and observations suggests that this consistency is driven by a realistic JRA55-do forcing. In general, in both hemispheres, an increase in resolution reduces SIE bias, with the exception of IAP-LICOM.

Figure 14 Modeled time evolution of annual mean Arctic and Antarctic (a, b) sea ice extent (10¹² m²) and (c, d) sea ice volume (10¹³ m³). The black lines are observations from the National Snow and Ice Data Center (NSIDC) and results the from Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) (Schweiger et al., 2011).

Sea ice volume (SIV) is a second independent measure of sea ice simulation performance that is representative of sea ice state, determined by the thermodynamical, optical, and dynamical properties of the ice itself, and the sensitivity of modeled SIV to model formulation is insightful for understanding the physical realism of the sea ice schemes. The modeled time evolution of annual mean Northern and Southern Hemisphere SIV is displayed in Fig. 14c and d. Relative to SIE, there is significant inter-model spread in SIV. With the exception of IAP-LICOM, all simulations exhibit a continuous decline in Northern Hemisphere sea ice volume during the satellite era, with trends similar to that of PIOMAS (Schweiger et al., 2011). Changes to SIV also do not appear to be consistent among the models when moving from low to high resolution: it increases for FSU-HYCOM and IAP-LICOM and decreases for NCAR-POP and AWI-FESOM. IAP-LICOM is a clear outlier in both hemispheres, which points to issues in the numerical implementation of its sea ice code at higher resolutions: sea ice thickness on average increases as the resolution is increased despite declining SIE. The IAP-LICOM sea ice model has only thermodynamic processes, and the absence of dynamic sea ice may lead to an excessive accumulation of volume.

3.6 Meridional overturning circulation

The Atlantic meridional overturning circulation (AMOC) is often quantified by a meridional overturning streamfunction with respect to depth, ψ_z at each latitude, defined as the meridional transport (Sv) across the basin above a constant depth z:

$${\mathit{\psi }}_z\left(yz\right)=\stackrel{\mathrm{‾}}{\underset{z}{\overset{\mathrm{0}}{\iint }}v\left(x^{\prime },y,z^{\prime },t\right)}\text{d}{z}^{\prime }\text{d}{x}^{\prime },$$

where v is the 4-D meridional velocity, the overbar denotes a time average (here annual mean), and the x integration covers the entire span of the Atlantic basin. The magnitude of the AMOC, or the AMOC transport, is then defined as the maximum of the streamfunction ψ_z with respect to z, representing the total northward transport above the overturning depth (approximately 1000 m in the Atlantic Ocean). This important measure of the AMOC is quantified and monitored by the RAPID moorings near 26.5^∘ N (e.g., Smeed et al., 2018).

The modeled annual mean AMOC transport across 26.5^∘ N from 1958 to 2018 is displayed in Fig. 15a, along with the RAPID results (2005–2017) for reference. The vertical distribution of the time mean zonal streamfunction in the Atlantic Ocean is shown later in Fig. 26. For the four models, the high-resolution simulations have a higher AMOC transport than their low-resolution counterparts: the mean AMOC transport in 2004–2018 is 7.8–14.9 Sv in the four low-resolution simulations versus 14.0–19.8 Sv in the four high-resolution simulations (17.2±1.5 Sv as observed by RAPID; e.g., Smeed et al., 2018). However, there is not an obvious difference in the AMOC sensitivity to forcing or in the trend between the low- and high-resolution experiments. The simulated time evolution of the AMOC transports at this latitude depends on the model, its parameterizations, resolution, and the number of spin-up cycles (Danabasoglu et al., 2016). Despite the fact that the simulations were only integrated over one JRA55-do cycle (1958–2018), all four low-resolution models show a similar multi-decadal variability, with a transport decrease from 1958 to the late 1970s, an increase from the late 1970s to late 1990s, and a decrease again thereafter. This multi-decadal variability is also present in the CORE-II simulations (Danabasoglu et al., 2016; Xu et al., 2019). In the high-resolution simulations, FSU-HYCOM, NCAR-POP, and AWI-FESOM exhibit a similar multi-decadal variability (of low transports in the late 1970s and high transports in the early 1990s) as seen in previous basin-scale simulations (e.g., Böning et al., 2006; Xu et al., 2013). The decline in the AMOC transport from the early 1990s to 2000s may be a consequence of the warming and reduced deep convection in the western subpolar North Atlantic that has been documented quite extensively (e.g., Häkkinen and Rhines, 2004; Yashayaev, 2007). The high-resolution AWI-FESOM has a strong AMOC decreasing trend (of 4–5 Sv) initially during the first 20 years and then has a time evolution of the AMOC that is very similar to that of FSU-HYCOM and NCAR-POP. The high-resolution IAP-LICOM shows an increasing AMOC transport from 1958 to the early 1990s (by 10 Sv) and a decrease thereafter (by 5 Sv), which is quite different from the other three models, but it does appear to capture the observed increase over the last few years.

Figure 15 Modeled annual mean Atlantic meridional overturning circulation (AMOC) transports at (a) 26.5^∘ N and (b) 34^∘ S; (c) global meridional overturning circulation (MOC) transports at 34^∘ S. The global MOC transport is negative because the northward-flowing Antarctic Bottom Water (AABW) is below the southward return flow (see Lumpkin and Speer (2007) and Talley (2013) for observed values). The solid black line in panel (a) is derived from the RAPID array measurements (Smeed et al., 2018).

AMOC transport in the South Atlantic Ocean has been quantified near 34^∘ S using several observational techniques including an expendable bathythermograph (XBT), Argo profiles, and moored current meter arrays in the western and eastern boundaries (Baringer and Garzoli, 2007; Dong et al., 2009, 2014, 2015; Garzoli et al., 2013; Goes et al., 2015; Meinen et al., 2013, 2018); these observations yield a time mean AMOC transport of about 14–20 Sv (see Xu et al., 2020a, for an in-depth discussion of circulation in the South Atlantic Ocean). The modeled temporal evolution of the AMOC transports at this latitude is overall similar to that at 26.5^∘ N in the North Atlantic (Fig. 15a, b). The range of the modeled time mean AMOC transport at 34^∘ S in the high-resolution simulations, from 14.7 Sv in FSU-HYCOM to 20.1 Sv in IAP-LICOM, is about the same as the observational range mentioned above.

The Antarctic Bottom Water (AABW) cell of the global overturning circulation (Talley, 2013) can be visualized as a streamfunction similar to Eq. (1), except that the zonal integration now spans across the full longitude circle. The transport associated with this cell is defined as the minimum of the streamfunction that is typically found near 3500 m (Lumpkin and Speer, 2007) and is a measure of the northward flow of the near-bottom water (AABW and/or Circumpolar Deep Water, CDW) from the Southern Ocean into the Atlantic, Pacific, and Indian Ocean that eventually upwells and returns southward at a shallower depth. Several observational estimates are provided based on hydrographic sections near 32^∘ S, but the uncertainty remains quite high (29±7.6 Sv in Talley, 2013). Figure 15c displays the modeled transport at 34^∘ S. The low-resolution simulations show a consistently lower transport, with the mean transport value for 2004–2018 ranging from 7.4 Sv in IAP-LICOM to 15.3 Sv in FSU-HYCOM. The low transport in IAP-LICOM is due to a long-term downward trend in this simulation: from 20 Sv in the early 1960 to 5 Sv in the late 2010s. The other three models, especially FSU-HYCOM, exhibit relatively stable transports in the last 30 years (of the integration). The high-resolution simulations show a much wider range, with time mean transports (2004–2018) ranging from 11.7 Sv in IAP-LICOM to 26.5 Sv in AWI-FESOM. The time evolution of the 34^∘ S transports in the high-resolution simulations typically shows an increase in the early stage of the integration followed by a gradual decrease afterward. But the timescale of the increase and decrease periods varies significantly between models. As a result, the transport is relatively stable in IAP-LICOM and AWI-FESOM for the last 30 years of the integration but continues to decrease slowly in NCAR-POP and FSU-HYCOM. Interestingly, despite the large difference in time mean transports, all four high-resolution models exhibit a similar interannual variability at this latitude for the last 30 years of the integration.

3.7 Drake Passage and Indonesian Throughflow transports

The Antarctic Circumpolar Current (ACC) is a strong oceanic current that flows eastward around Antarctica. It connects the Pacific, Atlantic, and Indian Ocean to its north and is the primary means of inter-basin exchange, enabling a truly global overturning circulation (e.g., Gordon, 1986; Schmitz, 1996; Talley, 2013). Accurate knowledge of ACC transport is fundamental to understand its influence on global circulation. Substantial observational efforts have been made toward quantifying ACC transport, especially in the Drake Passage, where the ACC is constricted between the Antarctica and the southern tip of South America. The major observations include (1) the International Southern Ocean Studies (ISOS) program in the late 1970s and early 1980s (Whitworth, 1983; Whitworth and Peterson, 1985), (2) the repeat hydrographic sections along the World Ocean Circulation Experiment (WOCE) line SR1b (e.g., Cunningham et al., 2003), (3) the repeat shipboard acoustic Doppler current profiler (ADCP) surveys that directly measure the current in the top 1000 m of the water column (Firing et al., 2011), (4) the DRAKE program that combines a moored current meter array (2006–2009) and satellite observations (e.g., Koenig et al., 2014), and (5) the cDrake program with a high-resolution bottom current meter mooring array and cPIES observations from 2007 to 2011 (e.g., Chidichimo et al., 2014; Donohue et al., 2016). The time mean, full-column transport estimates range from 134±15 Sv based on ISOS, to 141±2.7 Sv based on DRAKE program, and to 173.3±10.7 Sv based on the cDrake results.

Figure 16 Modeled annual mean volume transport of the ACC through the Drake Passage and of the Indonesian Throughflow (ITF) from the Pacific to Indian Ocean. The dotted lines in (a) represent the observational range between the canonical transport value (134 Sv) based on the ISOS program and 173 Sv based on cDrake observations. The dotted line in (b) is the observed ITF transport value (15 Sv) based on INSTANT observations.

The modeled annual mean Drake Passage transport is displayed in Fig. 16a together with the mean observational range (134–173 Sv). The transports in the low-resolution simulations differ from each other: IAP-LICOM and AWI-FESOM show a continuous decrease and FSU-HYCOM shows a continuous increase, whereas NCAR-POP shows an increase in the first 20 years, then levels off for the next 40 years of the integration. It is expected that the low-resolution models will be sensitive to different parameterizations and their ability to simulate eddy saturation (Gent and Danabasoglu, 2011; Munday et al. 2013). As a result, the time mean transport for the last 15 years of these low-resolution simulations is either near the edge or outside the above observational range: 132, 135, 173, and 190 Sv, respectively, for IAP-LICOM, AWI-FESOM, NCAR-POP, and FSU-HYCOM. The evolution of Drake Passage transport is slightly more similar in high-resolution simulations (compared to the low-resolution counterparts). Even at high resolution, mean ACC transport is sensitive to subgrid schemes (Pearson et al., 2017). Three models show a fast decrease in the first 10–15 years and then a very small decrease thereafter, except that FSU-HYCOM shows a small increase in the first 20 years and then levels off thereafter. The high-resolution models exhibit a higher interannual variability correlation than the low-resolution models. The ACC transports in NCAR-POP and AWI-FESOM of 120 and 130 Sv, respectively, are lower than the latest best mean transport estimate (Xu et al., 2020b) but still within the uncertainty range, while the transports in IAP-LICOM and FSU-HYCOM of 145 and 157 Sv, respectively, are within the mean observational range.

The long-term trend in modeled ACC transport merits some further discussions. Observations indicate that, as shown by a higher Southern Annular Mode index, the westerly winds in the Southern Ocean have been intensified since the 1970s. This increase in westerly winds is also present in the JRA55-do forcing, yet none of the high-resolution simulations exhibit a long-term upward trend (three models even have a slightly weakening trend). The ACC transport time series from the DRAKE program (Koenig et al., 2014) and from the WOCE line SR1b also show a stable ACC transport from early 1990s. Given that the ACC is wind-driven, a lack of dependence for ACC transport on zonal wind at a long-term timescale may at first be surprising, but it has also been shown that the strengthening of winds does generate more eddy activity along the path of the ACC, without necessarily changing its total transport (e.g., Hallberg and Gnanadesikan, 2006; Morrison and Hogg, 2013; Munday et al., 2013; Farneti et al., 2015). Following this argument, one may expect to see an increasing ACC transport in low-resolution, non-eddying simulations as long as the coefficient used in the eddy parameterization (i.e., GM) remains constant (Gent and Danabasoglu, 2011). Only FSU-HYCOM shows a long-term increase in the ACC transport, and two models (IAP-LICOM and AWI-FESOM) even show a long-term decrease. This may be a consequence of the eddy parameterization choice made by the different models, assuming that the ACC is quasi-equilibrated and that changes are not related to buoyancy changes in the Southern Ocean (e.g., rate of bottom water formation). There is, however, an initial condition transient in the buoyancy gradients across the ACC (Fig. 15c), and the MOC across the ACC is equilibrating for a significant fraction of the integration time: one cannot rule out the possibility that the ACC changes result from buoyancy drifts. All models have weaker Drake Passage transports at higher resolution, except for IAP-LICOM.

The Indonesian Throughflow (ITF) connects the tropical Pacific and Indian Ocean and provides a pathway for inter-basin exchange of heat and fresh water. Tropical Pacific and Indian Ocean water masses can go through the ITF and contribute to the upper AMOC limb in the South Atlantic Ocean as part of the global MOC via the Agulhas Current (e.g., Gordon, 1986; Schmitz, 1996; Talley, 2013). A good overview of the Indonesian Sea oceanography and the ITF observations can be found in Gordon (2005) and Gordon et al. (2010). The 3-year mean total ITF transport measured by the International Nusantara Stratification and Transport (INSTANT) program in 2004–2006 is 15±2.5 Sv (Sprintall et al., 2009; Gordon et al., 2010), about 30 % greater than the values of nonsimultaneous measurements made prior to 2000. ITF transport variability exhibits a close relationship to the phase of the El Niño–Southern Oscillation (Meyers, 1996).

The modeled annual mean ITF transport is displayed in Fig. 16b. The low-resolution simulations exhibit a wide range of ITF transports, ranging from about 5 Sv in the IAP model, to 11 Sv in the NCAR-POP model, and to 15 Sv in FSU-HYCOM and AWI-FESOM (30-year mean). ITF transports exhibit a gradual decrease in AWI-FESOM but are relatively stable in the other three models. In the high-resolution simulations, however, the models yielded a quite similar ITF transport, especially in the second half of the integration. All four models represented a similar mean transport of 15 Sv that is close to observations and all show a very similar variability. While large-scale forcing mechanisms are thought to set the basic level of ITF transport (Godfrey, 1989), the difference between the low-resolution and high-resolution simulations is likely impacted by the model ability to represent the small-scale bathymetry feature of the Indonesian Seas (Jochum et al., 2009); after all, the ITF enters the Indian Ocean through a number of narrow straits that are hard to represent in the low-resolution models. IAP-LICOM and NCAR-POP have low ITF transports at low resolution and appear to have stronger sensitivity to model resolution.

Figure 17 Top panels: mean 1993–2018 SSH and variance AVISO. Middle panels: difference between the mean modeled SSH and AVISO SSH. Bottom panels: modeled variance derived from 5 d average outputs. The low-resolution IAP-LICOM SSH variance was not provided.

4 Stationary ocean climate

4.1 Sea surface height (SSH) and eddy kinetic energy (EKE)

The impact on the circulation of increasing the horizontal resolution is twofold. First, the solution becomes more nonlinear and allows for a better representation of western boundary currents. Second, the first Rossby radius of deformation is resolved over most of the globe (Hallberg, 2013) and eddies are formed through barotropic and baroclinic instabilities, although higher vertical modes including submesoscale eddies are not often resolved (nor are they resolved by altimetry, i.e., AVISO) (see Stewart et al., 2017). Figure 17 displays the mean sea surface height bias with respect to AVISO and its 5 d variance over the 1993–2018 period. Overall, the large-scale patterns are well-represented globally, and one can clearly see the improved representation of the western boundary currents as resolution is increased (Gulf Stream and Kuroshio). In the North Atlantic, most coarse-resolution models to date have the tendency to exhibit an overshooting Gulf Stream and a poor representation of the North Atlantic Current (NAC) at the Northwest Corner. This is the case for three out of the four models. Instead of flowing northeastward along the continental rise past the Flemish Cap and continuing northwestward as in the observations (Rossby, 1996), the North Atlantic Current is strongly zonal in NCAR-POP, AWI-FESOM, and IAP-LICOM and does not turn northeastward at the Flemish Cap. This has been a long-standing issue for many ocean components of the CMIP climate models, and it does not necessarily improve as the computational mesh is refined. Increasing the horizontal resolution did improve the Gulf Stream separation (see Chassignet and Marshall, 2008, and Chassignet and Xu, 2017, for a review) in all models but not necessarily the representation of the Northwest Corner circulation (Fig. 17; see Fig. 1 in Chassignet et al., 2020a, for details). FSU-HYCOM is the only model that has a good representation of the Northwest Corner in both the low- and high-resolution experiments. Since all the models use the same atmospheric forcing dataset, the difference is solely due to the numerical and physical choices made by each modeling group.

As expected, there is a significant increase in the SSH variance as resolution is increased, and the eddying solution SSH variance maps are much closer to the observations (top panels) than their low-resolution counterparts. In the high-resolution experiments (∼ 0.1^∘), the variability is, however, still lower than observed, especially in the gyre interiors and in the experiments that use relative winds (NCAR-POP, AWI-FESOM, and IAP-LICOM). This underestimation is thus partly a consequence of the well-known eddy-killing effect, which results from taking into account the shear between atmospheric wind and ocean current when computing the wind stress and which can reduce the kinetic energy by as much as 30 % (see discussion in Sect. 3.1). However, the use of absolute wind in FSU-HYCOM is not sufficient to raise the level of surface variability to that of the observations, and Chassignet and Xu (2017) argue that one actually needs to significantly increase the resolution (∼ 0.01^∘) in order to resolve the submesoscale instabilities that can energize the mesoscale (Callies et al., 2016) and therefore enhance eddy kinetic energy comparable to the mesoscale AVISO observations. It is more physical to take into account the vertical shear between atmospheric winds and ocean currents when computing the wind stress (see Renault et al., 2020, for a review) as it allows for a better representation of western boundary current systems (Ma et al., 2016), especially the Agulhas Current retroflection and associated eddies (Renault et al., 2017). In FSU-HYCOM, which uses absolute wind, the Agulhas eddies are too regular and follow the same pathway. The use of relative winds not only increases the eddy decay, but also impacts the location of the Agulhas retroflection and where the eddies are formed. The three simulations with relative winds have a better representation of the Agulhas eddy pathways.

Figure 18 Surface EKE from (a) AVISO and (b) drifter observations. Modeled surface EKE calculated from 5 d average fields for the high-resolution (c) FSU-HYCOM, (d) NCAR-POP, (e) AWI-FESOM, and (f) IAP-LICOM.

The eddy kinetic energy maps for the high-resolution simulations are displayed in Figs. 18–21 for the surface, at 700 m, and at 1000 m. Observed surface EKE maps are either derived from altimetry (Fig. 18a) or drifters (Fig. 18b). Because of the irregular sampling, this necessarily implies some type of smoothing in space and time. In the case of the AVISO along-track altimetry observations, they are optimally interpolated on a 0.25^∘ grid, which filters scales less than 150 km (due to track separation and measurement noise and errors) and timescales less than 10 d (repeat cycle of the altimeters). There is a significant reduction in EKE when computing it using 5 or 10 d average outputs (as opposed to online at every time step) as shown in Fig. 19 for FSU-HYCOM and NCAR-POP. In addition to a significant EKE reduction in the most active regions, the time averaging removes much of the small-scale variability associated with inertial motions and ageostrophic effects. There is not a big difference between the 5 d and the 10 d maps, and for the rest of this section we will use 5 d average model outputs when comparing the EKE to observations. Overall, as for the SSH variability, the EKE is larger in FSU-HYCOM because absolute winds are used to force the models. As already discussed, the use of relative winds does improve the pathway for the Agulhas Current eddies, but in Fig. 18 there is also a significant reduction of EKE in the ACC, and it also suppresses variability in many areas. This is especially true for the Indian Ocean, in the tropics, and west of the Hawaiian Islands.

Most model–observation comparisons usually focus on the surface fields because of very sparse spatiotemporal sampling at depth (see, for example, the EKE derived at 700 m from several years of SOFAR float measurements in the Gulf Stream region, Richardson, 1993, or the EKE distribution from ARGO floats at 1000 m, Ollitrault and Colin de Verdière, 2014). Scott et al. (2010) evaluated the total kinetic energy simulated in eddying ocean models (on the order of 0.1^∘ horizontal resolution) relative to moored current records. They found that the models agreed within a factor of 2 above 3500 m of depth and within a factor of 3 below 3500 m of depth. Thoppil et al. (2011) show that the surface and abyssal EKE increases with resolution and that better upper-ocean EKE allows strong eddy-driven abyssal circulation. This also means that, while the EKE at depth in the high-resolution experiments is a significant improvement over the quasi-nonexistent EKE of the low-resolution simulations, it is still significantly less than very limited observations (Richardson, 1993; Ollitrault and Colin de Verdière, 2014). Significantly higher resolution may be necessary in order to obtain a level of EKE close to the observations (see Chassignet and Xu, 2017, for a discussion). When using relative winds at this resolution as in NCAR-POP, AWI-FESOM, or IAP-LICOM, there is very little EKE at depth when compared to FSU-HYCOM forced with absolute winds (Figs. 20 and 21).

Figure 19 Comparison of surface EKE calculated from the total integrated kinetic energy, 5 d average fields, and 10 d average fields in the high-resolution FSU-HYCOM and NCAR-POP.

4.2 Sea surface temperature (SST) and sea surface salinity (SSS) deseasonalized variance

Surface temperature and salinity show sensitivities to resolution that are largely consistent with those described above for SSH in terms of both mean state and variability. All models show greatly enhanced SST variance (computed from deseasonalized monthly values spanning 1980–2018) at high resolution in regions of high mesoscale eddy activity, including the Southern Ocean, the Agulhas Current retroflection, the Brazil–Malvinas Confluence region, the Kuroshio and its extension across the North Pacific, and the Gulf Stream and its NAC extension (Fig. 22). A global, satellite-based SST dataset spanning 1982–2018 provides some measure of ground truth for comparison, but the sampling limitations of microwave and infrared measurements of SST restrict the observed estimate to a (1∕4^∘) grid that is 2–3 times coarser than the high-resolution model grids considered here (OISST.v2; Reynolds et al., 2007; Banzon et al., 2014).

Figure 20 EKE at 700 m calculated from 5 d average fields in (a) FSU-HYCOM, (b) NCAR-POP, (c) AWI-FESOM, and (d) IAP-LICOM.

Figure 21 EKE at 1000 m calculated from 5 d average fields in (a) FSU-HYCOM, (b) NCAR-POP, (c) AWI-FESOM, and (d) IAP-LICOM.

Figure 22 Modeled and observed (OISST) deseasonalized SST variance.

As was noted for SSH, the high-resolution FSU-HYCOM simulation generally exhibits the highest SST variance in the midlatitude gyre interiors compared to other models, but it has lower SST variance in strongly eddying regions such as the Agulhas retroflection, Kuroshio, and NAC. Overall, the high-resolution FSU-HYCOM shows the closest match to the observed benchmark. The improved structure of SST variance in the subpolar Atlantic in FSU-HYCOM is related to the improved NAC pathway in that model (discussed above). Compared to OISST.v2, the low-resolution models generally underestimate SST variance, while the high-resolution models tend to overestimate the variance, particularly in the western boundary currents, Southern Ocean, and Agulhas retroflection region. The resolution sensitivity is most dramatic in the NCAR-POP model, whose low-resolution version shows the weakest SST variance of any of the simulations considered here. The overestimation of SST variance has been noted before in eddy-resolving coupled simulations (e.g., Small et al., 2014), but it is not clear whether this discrepancy is attributable to deficiencies in the models or in the observation-based estimates, which cannot fully resolve ocean mesoscale variability. The absence of a dynamic (and ocean mesoscale-resolving) atmosphere in these simulations may partially explain the overly high SST variance insofar as important eddy-damping processes (Ma et al., 2016) are missing in this experimental framework. It is interesting that all high-resolution simulations generate more SST variance in the Kuroshio region than seen in OISST.v2, and yet they all show less variance than observed in the northeast Pacific and along the southeastern branch of the North Pacific gyre (near Hawaii). There is an indication of slightly enhanced ENSO-related tropical Pacific SST variance when going from low to high resolution, but there is not a systematic change in the spatial structure of this variance, which appears to depend mainly on model formulation. The representation of high extratropical SST variance along the eastern boundaries of each of the major ocean basins shows robust improvement with resolution across the different models. The high-resolution versions all show higher (and more realistic) variance in the following locations: the Benguela Current region in the eastern South Atlantic; the coastal region off the west coast of North America and Baja; the Leeuwin Current region off the west coast of Australia; and the Canary Current region off the west coast of North Africa. Improvements in simulated SST variability in these regions are more apparent than improvements in the mean state (Fig. 10). The overall picture of greatly enhanced SST variance at eddying resolution is an important, but not unexpected, result that has significant implications for climate modeling because resolving air–sea interactions at the ocean mesoscale has been shown to result in qualitatively different coupled model behavior (e.g., Bryan et al., 2010; Ma et al., 2016, 2017; Laurindo et al., 2019).

Figure 23 Modeled and observed deseasonalized SSS variance.

While surface salinity generally exhibits globally enhanced variance at high resolution, particularly in the eddy-rich locations already highlighted above, the SSS variance plots are striking in that they show considerably more sensitivity to model formulation than to model horizontal resolution (Fig. 23). The increase in SSS variance at both resolutions across model systems – from FSU-HYCOM, to NCAR-POP, to AWI-FESOM, and to IAP-LICOM – is likely related to the steady increase in the surface-salinity-restoring timescale (Table 1). In the absence of reliable long-term global SSS observations, it is difficult to say which model formulation is best or whether high resolution improves the simulation of surface salinity variability.

4.3 Upper-ocean heat and salinity content mean biases and deseasonalized variance

The temperature and salinity distributions in the upper 700 m are less constrained to observations than the SST (which is restored to time-varying observed values through the sensible heat flux) and the SSS (which is restored to observed climatology through the artificial salinity-restoring flux). Therefore, vertically averaged upper-ocean heat content (UOHC) and salt content (UOSC) biases tend to be larger than surface biases, particularly at low latitudes (see Figs. 10, 24, and 25). The UK Met Office EN4 ocean analysis (Good et al., 2013) provides an observational benchmark for UOHC and UOSC mean and variance, but again considerable caution is warranted when interpreting model–observation discrepancies. The spatial resolution of this analysis is only 1^∘, and it relaxes to climatology in data-sparse regions, which are extensive in the pre-Argo era (prior to 2000).

Both FSU-HYCOM and NCAR-POP show a significant and nearly ubiquitous mean bias reduction for UOHC and UOSC when resolution is increased (Figs. 24 and 25). In contrast, AWI-FESOM and IAP-LICOM show mixed results, with bias reduction in some regions (e.g., the Southern Ocean) offset by bias increase elsewhere. Regions of high observed UOHC variance include all the regions of high observed SST variance mentioned above as well as some variance hot spots in the subtropical western Pacific, the Tasman Sea, the subtropical south Indian Ocean, and the western tropical Atlantic (Fig. 24). The high-resolution models tend to show improved representation of these subsurface variance hot spots and enhanced variance in high EKE regions (western boundary currents, Agulhas, Southern Ocean, etc.). The high-resolution IAP-LICOM stands out for its globally high (perhaps unrealistically high) UOHC variance. This upper-ocean variance bias in IAP-LICOM is even more evident in UOSC (Fig. 25), although the limited sampling of 0–700 m salinity in the EN4 product implies that the observed variance estimate is almost certainly a gross underestimate. There is a hint (especially apparent in NCAR-POP and IAP-LICOM) of higher resolution resulting in enhanced salinity variance along subtropical cell spiciness pathways from the eastern extratropical Pacific to the tropical western Pacific (e.g., Yeager and Large, 2004).

Figure 24 Top panels: mean 1980–2018 0–700 m temperature and variance from EN4. Middle panels: difference between the mean modeled 1980–2018 0–700 m temperature and EN4. Bottom panels: deseasonalized 0–700 m temperature modeled variance.

Figure 25 Top panels: mean 1980–2018 0–700 m salinity and variance from EN4. Middle panels: difference between the mean modeled 1980–2018 0–700 m salinity and EN4. Bottom panels: deseasonalized 0–700 m salinity modeled variance.

4.4 Zonal AMOC mean and variance

Although the magnitude and temporal evolution of AMOC transport, as shown in Fig. 15a and b, differ significantly between different models and different horizontal resolutions, their time mean spatial structure and variance are quite similar among all the simulations as shown in Fig. 26 for the 1980 to 2018 period. All simulations show a positive upper cell and a negative lower cell. In the low-resolution simulations, the upper cell exhibits maximum transport near 1000 m and 35–40^∘ N. The lower cell has its maximum transport near 3500–4000 m and is much weaker overall (∼ 2 Sv), except for IAP-LICOM (∼ 6 Sv in the subtropical North Atlantic). In high-resolution simulations, the upper cell extends deeper in all four models, which is more consistent with the observations, and the lower cell has a similar magnitude (2–4 Sv). This can be seen more clearly in Fig. 27, which compares the model AMOC profiles for 2004–2018 to the RAPID results near 26.5^∘ N. South of about 20^∘ N, the high-resolution AWI-FESOM exhibits a weak positive cell near the bottom that is unrealistic and is different from other models.

Figure 26 Top panels: modeled time mean (1980–2018) Atlantic meridional overturning streamfunction ψ_z(yz) in the four models (low and high resolution). Lower panels: modeled standard deviation (1980–2018) of the Atlantic overturning streamfunction ψ_z(yz) derived from monthly mean fields.

The AMOC standard deviation shows a similar pattern between different resolutions and different models (Fig. 26). The variability is overall highest near the Equator (IAP-LICOM has the strongest variability near 10^∘ N that differs from the other three). This is due to the large seasonal variability of the AMOC associated with the migration of the Intertropical Convergence Zone (ITCZ) and the changes in wind patterns and Ekman transport (Xu et al., 2014). The standard deviations based on annual means do not show such a maximum near the Equator (see Fig. 7 in Hirschi et al., 2020, for results based on different averaging scales). The variability is much weaker beyond the equatorial region. In the low-resolution simulations, the standard deviation is typically 2 Sv and has a slight elevation near 40^∘ N. In the high-resolution simulations, the standard deviation is about 3 Sv and clearly shows a secondary maximum near 40^∘ N. The difference between the low- and high-resolution simulations near 40^∘ N highlights the impact of the meandering and mesoscale eddy variability of the Gulf Stream on AMOC variability.

Northward heat transport is usually correlated with higher AMOC transport in the Atlantic Ocean (e.g., Msadek et al., 2013; Xu et al., 2016). In the high-resolution experiments, there is also an overall reduction in subtropical Atlantic upper-ocean cold bias (Fig. 24), which can have an impact as shown by Msadek et al. (2013). Thus, it is not surprising that the high-resolution simulations have a higher heat transport in the Atlantic than the low-resolution counterparts (Fig. 28). The maximum northward heat transport at 20–30^∘ N is about 1.0 PW in the high-resolution models, 0.6–0.8 PW in low-resolution models, and 1.25–1.3 PW in observations (Johns et al., 2011; Trenberth et al., 2019). In other basins such as the Indo-Pacific and the Southern Ocean, higher resolution does not necessarily lead to higher heat transport (Fig. 28a–c). But in general, the meridional heat transport in high resolution is closer to observations than the low-resolution simulations (Fig. 28).

Figure 27 AMOC profile at 26^∘ N.

Figure 28 Modeled time mean meridional heat transport in PW over 1980–2018: (a) global, (b) Atlantic Ocean, and (c) Indo-Pacific Ocean. Black lines are the latest observational results from Trenberth et al. (2019) with uncertainties.

4.5 Northern and Southern Hemisphere sea ice mean and variance

Figures 29–32 show 1980–2018 March and September modeled mean sea ice concentration (SIC) and sea ice thickness (SIT) together with the mean 1980–2018 SIC from passive microwave satellites and the mean 2003–2007 SIT from ICESat. As expected, given the relationship between sea ice extent and atmospheric temperature, the observed spatial pattern of winter SIC (averaged from 1980 to 2018) in the Northern Hemisphere (Fig. 29) is well-simulated in all the models at both low and high resolutions though with greater consistency among the high-resolution models. The position of the marginal ice zone (MIZ) – or areas where the SIC is between 15 % and 80 % – varies between models and resolutions. With the SIC close to 100 % in the interior Arctic in winter, the variance in modeled SIC (Figs. S1–S8) is small in March. The general pattern of the March Northern Hemisphere SIT is consistent between models, save for the IAP-LICOM, with increasing SIT towards the Canadian Archipelago. At both resolutions, AWI-FESOM presents slightly thicker sea ice throughout most of the Arctic Basin. The high-resolution IAP-LICOM is unrealistic, with very high sea ice thicknesses surrounding the pole. Figure 30 repeats Fig. 29, but for September. The comparison is similar: there are large inter-model differences in the MIZ, but there is a generally consistent pattern of increasing ice concentrations and thickness from the pole toward the Canadian Archipelago and northern Greenland. The location of the summer MIZ in FSU-HYCOM and NCAR-POP is too far away from the coast when compared to the observations and AWI-FESOM. Overall, there is a broad inter-model disagreement in SIT. The SIT variance (Figs. S1–S8) is higher in summer, especially in the shelf regions in experiments at both resolutions. While there is some sensitivity to model resolution for the simulated sea ice thickness, the difference due to changes in the resolution is smaller than the difference between the model results and the observations. Since the SIT observations shown in Figs. 29–32 are 2003–2007 averages, one may question the value of comparing the 2003–2007 ICESat observations to 1980–2018 averaged modeled fields. For a more direct comparison, the SIT modeled fields averaged from 2003 to 2007 are shown in Figs. S9 and S10 together with the ICESat observations when available, but they are almost identical to the 1980–2018 SIT fields shown in Figs. 29–32.

Figures 31–32 reproduce Figs. 29–30 for the Southern Hemisphere. The geographically simpler observed SIC in winter is well-represented by the models, with the exception of FSU-HYCOM, which simulates large regions of low ice concentration within the Antarctic pack. There is a clear difference in the modeled MIZ across the simulations: for example, it is wider in the low-resolution IAP-LICOM experiments when compared to the high-resolution experiment and vice versa for NCAR-POP. In the AWI-FESOM simulation, the MIZ is much narrower. The MIZ width is an important observable for explaining Southern Hemisphere sea ice concentration and thickness, largely determined by the presence of ocean surface waves (Kohout et al., 2014; Horvat et al., 2020). As none of the models allow for the propagation of ocean waves into the ice, the feedback between wave propagation and sea ice cover could not take place. In summer, SIT and SIC are poorly represented in FSU-HYCOM, NCAR-POP, and IAP-LICOM, which, at high resolution, produce an MIZ that is too wide (FSU-HYCOM, NCAR-POP) or no MIZ (IAP-LICOM) and inaccurately low representations of SIT, with most regions thinner than 20 cm. In general, observations indicate that areas of the most compact sea ice (Weddell Sea) are not those of the thickest ice (western Antarctic sector) in summer. None of the simulations adequately represent the observed relatively thick and noncompact sea ice in the western Antarctic. The high summer sea ice concentration in the Weddell and Ross Sea is only reasonably simulated in AWI-FESOM but with ice that is too thick. Although there is a high sensitivity to changing resolution in the Southern Hemisphere, the low- and high-resolution model pairs tend to agree more than the ensemble of low- or high-resolution models, indicating that sea ice model formulation differences may exceed the effects of resolution in these patterns. The models also show some sensitivity to the model resolution for the variance of the sea ice concentration and thickness in the Southern Hemisphere, especially in FSU-HYCOM. However, the difference between models is considerably larger. As IAP-LICOM significantly overestimates the sea ice thickness in its high-resolution setup, its variance is also much higher than the other models (Figs. S1–S8).

Figure 29 Northern Hemisphere winter mean (March) sea ice concentration (SIC) and thickness (SIT). Top panels: 1980–2018 SIC from passive microwave satellites and 2003–2007 SIT from ICESat. Lower panels: 1980–2018 modeled SIC and SIT.

Figure 30 Northern Hemisphere summer mean (September) sea ice concentration (SIC) and thickness (SIT). Top panel: 1980–2018 SIC from passive microwave satellites (summer SIT not available from ICESat). Lower panels: 1980–2018 modeled SIC and SIT.

Figure 31 Southern Hemisphere winter mean (September) sea ice concentration (SIC) and thickness (SIT). Top panels: 1980–2018 SIC from passive microwave satellites and 2003–2007 SIT from ICESat. Lower panels: 1980–2018 modeled SIC and SIT.

Figure 32 Southern Hemisphere summer mean (March) sea ice concentration (SIC) and thickness (SIT). Top panels: 1980–2018 SIC from passive microwave satellites and 2003–2007 SIT from ICESat. Lower panels: 1980–2018 modeled SIC and SIT.

5 Summary and discussion

The overall goal of this paper is to assess the robustness of climate-relevant improvements in ocean simulations (mean and variability) associated with moving from coarse (∼ 1^∘) to eddy-resolving (∼ 0.1^∘) horizontal resolutions. It also lays out a set of basic large-scale diagnostics for assessing the relative quality, variability, and sensitivity of high-resolution versus low-resolution ocean and sea ice models. The emphasis is on the key metrics used in climate modeling – SST, OHC, sea level, salinity, sea ice extent and volume, and circulations that tend to have global impacts (MOC, ACC, ITF) on the modeled climate. Here these metrics are assessed in a suite of four pairs of low-resolution–high-resolution ocean and sea ice models forced with the latest JRA55-do dataset (Tsujino et al., 2018). These results will provide a useful baseline for future process-focused analyses and ocean model development activities at diverse resolutions.

On the whole, the biases in the low-resolution variants are familiar and are similar in these models to those found when assessing sensitivity to forcing products (Tsujino et al., 2020). Gross features of the bias patterns in low-resolution models – position, strength, and variability of western boundary currents, equatorial currents, and ACC – are significantly improved in the high-resolution models. However, despite the fact that the high-resolution models “resolve” these features, the improvements in temperature and salinity are inconsistent among the different model families, and some regions show increased bias over their low-resolution counterparts. SSH variability and near-surface EKE are significantly – even qualitatively – improved in all high-resolution models over their low-resolution counterparts, although all of these models still underpredict the observed SSH variability and EKE, particularly in the ocean interior, which indicates a need for further refinements in resolution (Chassignet and Xu, 2017), and improvements in less dissipative subgrid schemes for high-resolution models (Pearson et al., 2017) are needed. The results in coupled models in the HighResMIP ensemble (Haarsma et al., 2016) show similar improvements in SSH and SST variability as well as EKE. Considerable differences in the high-resolution models used here were associated with the use of relative winds versus absolute winds. Renault et al. (2020) show that using relative winds contributes a feedback that tends to reduce EKE, and the models here are consistent with a 30 % or so reduction of EKE when relative winds are used.

One interesting aspect of the high-resolution models versus the low-resolution models was that the interannual variability in ACC and ITF transport as well as the AMOC was more consistent among the high-resolution models than among the low-resolution models. ITF transport was especially in agreement among the high-resolution models, indicating that better representation of the passageways through the Indonesian archipelago is critical. Consistency in all of these transports potentially indicates that higher-resolution models are needed to represent process variability, which may explain some of the past difficulties in comparing the magnitude of these phenomena across coarse-resolution models. However, the mean ACC transport and MOC strength were not in greater agreement among the high-resolution than the low-resolution models, which means that more work remains in evaluating sensitivity to numerics and subgrid-scale schemes for high-resolution models. Furthermore, Danabasoglu et al. (2016) note that low-resolution models come into greater agreement in AMOC variability after more cycles of the CORE forcing – this comparison was limited by the cost of the high-resolution models to only a single cycling of the forcing. It is unclear if the high-resolution models are in greater agreement only in the first cycle or generally. The short duration of a single forcing cycle limits the comparison of the decadal changes that are emphasized in Danabasoglu et al. (2016), so the improved agreement discussed here among high-resolution models is year by year rather than decade by decade. Nonetheless, the high-resolution models had systematically stronger and more variable AMOC, which was in better agreement with observations in both the maximum overturning and profile than the low-resolution models.

From a climate modeling perspective, ocean heat content, sea level, and sea ice stability are key metrics. There is some indication that the high-resolution models may warm more quickly below 700 m than the low-resolution models, indicating errors in parameterizations of vertical eddy heat transport. Griffies et al. (2016) found similar sensitivity to resolution in a model hierarchy, which they explained as resolved eddies versus parameterized eddies affecting vertical heat transport. However, warming between 0 and 700 m and global mean sea level rise were not systematically different across model families between the low-resolution and high-resolution models. Additional insights could be gained by computing vertical heat and salt budgets as in Griffies et al. (2015) and von Storch et al. (2016). Regional sea level rise, however, was significantly more variable in the high-resolution models wherein mesoscale eddies and the variability of western boundary and equatorial currents make an impact on regional sea level rise. Thus, it is critical to use high-resolution models to assess the ocean dynamic sea level (Gregory et al., 2019) contribution to regional sea level rise. Generally, there is a highly variable representation of sea ice with the exception of Arctic sea ice concentration – with high inter-model differences in sea ice thickness and volume, as well as generally inconsistent representation of sea ice in the Southern Hemisphere and sea ice thickness in the Arctic. Inter-model differences are larger than improvements with changing resolution – pointing to the known role of atmospheric forcing in setting modeled sea ice states and to the poor representation of sea ice physics and associated parameterizations. Sea ice concentration bias reduction is achievable, for example, through the inclusion of lateral melting and feedbacks between ice fracture and sea ice melting (Horvat and Tziperman, 2015; Roach et al., 2018, Bateson et al., 2020), but the wide spread in sea ice thickness and poor representation of Antarctic sea ice are known in CMIP-class models (Roach et al., 2020; Notz et al., 2020, Shu et al., 2020) and indicate that a fundamental overhaul of sea ice modeling may be needed.

Low-resolution versus high-resolution comparisons are often motivated by identifying persistent patterns across multiple models of bias or improvements with resolution. However, as is apparent in many metrics, there is no overall consistency among either the low-resolution or high-resolution ensemble from which to draw simplistic conclusions. This paper therefore does not dwell on the detailed differences among these models' numerics and parameterizations but is instead intended to serve as a benchmark for future studies comparing and improving different schemes in any of these models or similar ones. The numerics and parameterizations of these models are in continual development in both low-resolution and high-resolution versions. Ideally, however, one should strive to minimize the number of changes as the resolution is increased to ensure traceability. Other modeling centers have expressed interest in participating in this protocol in the future, and this paper provides a basis for comparison.

Finally, the models used in the present study are not coupled to active atmospheres, biospheres, land models, or land ice models. Thus, many important feedbacks and aspects of climate modeling are not addressed here. However, these same feedbacks complicate diagnosis of the source of model biases and resolution sensitivity, which is where this study intends to contribute.