The very-high-resolution configuration of the EC-Earth global model for HighResMIP

We present the very-high-resolution (VHR) version of the EC-Earth global climate model, EC-Earth3P-VHR, developed for HighResMIP. The model features an atmospheric resolution of ∼16 km and an oceanic resolution of $\mathrm{1}/\mathrm{12}$° (∼8 km), which makes it one of the finest combined resolutions ever used to complete historical and scenario-like CMIP6 simulations. To evaluate the influence of numerical resolution on the simulated climate, EC-Earth3P-VHR is compared with two configurations of the same model at lower resolution: the ∼100 km grid EC-Earth3P-LR (LR) and the ∼25 km grid EC-Earth3P-HR (HR). Of the three configurations, VHR shows the smallest drift in the global mean ocean temperature and salinity at the end of a 100-year 1950s control simulation, which points to a faster equilibrating phase than in LR and HR. In terms of model biases, we compare the historical simulations against observations over the period 1980–2014. In contrast to LR and HR, VHR shows a reduced equatorial Pacific cold tongue bias, an improved Gulf Stream representation with a reduced coastal warm bias and a reduced subpolar North Atlantic cold bias, and more realistic orographic precipitation over mountain ranges. By contrast, VHR shows a larger warm bias and overly low sea ice extent over the Southern Ocean. Such biases in surface temperature have an impact on the atmospheric circulation aloft, connected with a more realistic storm track over the North Atlantic yet a less realistic storm track over the Southern Ocean compared to the lower-resolution model versions. Other biases persist or worsen with increased resolution from LR to VHR, such as the warm bias over the tropical upwelling region and the associated cloud cover underestimation, a precipitation excess over the tropical South Atlantic and North Pacific, and overly thick sea ice and an excess in oceanic mixing in the Arctic. VHR shows improved air–sea coupling over the tropical region, although it tends to overestimate the oceanic influence on the atmospheric variability at midlatitudes compared to observations and LR and HR. Together, these results highlight the potential for improved simulated climate in key regions, such as the Gulf Stream and the Equator, when the atmospheric and oceanic resolutions are finer than 25 km in both the ocean and atmosphere. Thanks to its unprecedented resolution, EC-Earth3P-VHR offers a new opportunity to study climate variability and change of such areas on regional and local spatial scales, in line with regional climate models.

1 Introduction

Interest in high-resolution modeling has soared in the past years, especially thanks to large European research projects and initiatives such as PRIMAVERA (https://www.primavera-h2020.eu/, last access: 20 June 2024, PRIMAVERA and the European Commission, 2015), nextGEMS (https://nextgems-h2020.eu/, last access: 20 June 2024, Hohenegger et al., 2023, Rackow et al., 2024), EERIE (https://eerie-project.eu/, last access: 20 June 2024), and Destination Earth (https://destination-earth.eu/, last access: 20 June 2024, Hoffmann et al., 2023). Broadly, these projects seek to build the next generation of high-resolution global climate (or Earth system) models capable of representing climate phenomena with unprecedented accuracy to simulate and predict regional climate, guide policymaking, and provide relevant climate information to end users. Thanks to these efforts, high-resolution models at resolutions of 25–50 km or even finer have been proven to lead to reduced biases in the simulated climate (see Introduction in Moreno-Chamarro et al., 2022, for a review) and to a better representation of, for example, tropical cyclones (Roberts et al., 2020a; Vidale et al., 2021; Zhang et al., 2021), storm tracks (e.g., Hodges et al., 2011), the Intertropical Convergence Zone (ITCZ; e.g., Doi et al., 2012; Tian and Dong, 2020), and the Gulf Stream and associated air–sea interactions (e.g., Kirtman et al., 2012; Bellucci et al., 2021) compared to standard-resolution models (hereafter, ∼100 km grid). An extensive review of the benefit of high-resolution modeling can be found in Haarsma et al. (2016), Hewitt et al. (2017), Roberts et al. (2018b), and Czaja et al. (2019). However, increased model resolution alone is not always the answer: for example, persistent, well-known biases in clouds and radiation can be insensitive to an increase in atmospheric resolution from a ∼100 km grid to a 25–50 km grid (Moreno-Chamarro et al., 2022). Inadequate model physics or insufficient tuning can thus mask or negate the benefits of increased resolution.

High-resolution modeling faces additional challenges. One is the high computational cost of running the simulations, and another related one is the difficulty of achieving high throughput due to the loss of efficiency with increasing parallelization. These issues have gradually improved thanks to steady increases in supercomputing power and parallel enhancements in model efficiency to leverage that power. The community trusts high-performance computing (HPC) to increase the performance of climate models, developing different approaches to speed models up. These approaches can go from improving the traditional parallelization algorithms (Tintó et al., 2019) or reducing the accuracy of the variables from double to single precision (Váňa et al., 2017; Tintó Prims et al., 2019) to increasing the input/output throughput of complex model configurations (Yepes-Arbós et al., 2022; Sarmany et al., 2024). Faster models are also needed to complete, in a reasonable time, the tuning and the spin-up phases, which for a high-resolution model can be extremely costly. The demand for high efficiency in high-resolution modeling has therefore accelerated the development and implementation of new modeling strategies to ensure optimal use of computing resources.

High-resolution models also need to find a fair compromise between the resolutions of the different climate components, which, sometimes, can be very disparate – for example, an eddy-rich ocean model (∼10 km grid) coupled to a 25 km, 50 km, or even coarser-grid atmosphere model (e.g., Gutjahr et al., 2019; Rackow et al., 2019; Semmler et al., 2020). Tsartsali et al. (2022), for example, reported increased ocean–atmosphere coupling strength and better agreement with reanalysis and observations over the Gulf Stream when both the ocean and atmosphere resolutions were increased to a comparable ∼25 km grid at least. Moreton et al. (2021) showed a degraded representation of the air–sea interaction at increased oceanic resolution but a constant atmospheric resolution. Similarly, Ma et al. (2017) found that the mesoscale ocean temperature affects the storm track over the Pacific only when the atmospheric model resolution is enough to resolve the small-scale diabatic heating. Finally, Rai et al. (2023) described disproportionate eddy killing when a coarse 200 km wind forcing was used to force a finer (∼10–25 km) ocean compared to the case with similar grid sizes. These results of these studies thus advocate for a similar resolution in both the atmosphere and ocean.

Sometimes, high-resolution modeling relies on a single model component, either atmospheric-only (Baker et al., 2019) or ocean-only configurations (e.g., Biastoch et al., 2021), or on regional models (e.g., Woollings et al., 2010; Ma et al., 2017) as in CORDEX (Jacob et al., 2014) for hypothesis testing and downscaling climate projections. Such configurations, however, lack global energy constraints, remote influences, and, potentially, key feedbacks rectifying the mean state. These models are also limited by the boundary conditions, which often are derived from coarser (∼100 km) global models and can present biases in their mean climate that might be absent or much reduced at a higher resolution; these biases might then be passed onto the single-model configurations. For example, an overly smooth Gulf Stream temperature gradient, an incorrect separation, or the lack of mesoscale representation in ocean temperatures can impact the response of the atmospheric circulation aloft (e.g., Ma et al., 2017; Lee et al., 2018). Low-resolution and high-resolution global models can also respond differently to climate change: for example, the northward shift and strong surface warming of the Gulf Stream projected by the eddy-rich configuration of the HadGEM3-GC3.1 model for the 21st century are absent in the lower-resolution model versions (Moreno-Chamarro et al., 2021). Associated with this, the increase in winter precipitation is similarly much larger over Europe at the highest resolution than at any lower one, which reinforces the idea that the response of the atmosphere is strongly sensitive to the boundary conditions. These findings put a limit on our confidence in single-model configurations and regional models, since they lack a global dynamical response.

As a response to the listed challenges, we present the eddy-rich version of the EC-Earth climate model for PRIMAVERA and HighResMIP. This is likely one of the finest combined horizontal resolution global models ever used to complete CMIP-like simulations, with a nominal resolution of about 10–15 km; it also has the additional advantage that the resolution is comparable in both the atmosphere and ocean–sea ice, which allows the atmosphere to “see” the fine-scale forcing from the ocean with minimal information lost from interpolation. In this paper, we describe the model configuration and the developments in model efficiency (Sect. 2), as well as the main characteristics of its climate for the period 1980–2014 compared to observations (Sect. 3).

2 Model description and experimental setup

2.1 Model description

All HighResMIP contributions with the EC-Earth global coupled climate model have been performed with its version 3.2.2, developed within the PRIMAVERA project (EC-Earth3P). The model consists of atmosphere, ocean, and sea ice components. The atmosphere model is based on the ECMWF Integrated Forecasting System (IFS) in the 36r4 cycle (based on IFS system 4, https://www.ecmwf.int/sites/default/files/elibrary/2011/11209-new-ecmwf-seasonal-forecast-system-system-4.pdf, last access: 8 November 2024). A detailed account of the changes introduced in this cycle can be found on the ECMWF website (https://confluence.ecmwf.int/display/FCST/Implementation+of+IFS+Cycle+36r4, last access: 20 June 2024). The very-high-resolution version of the model, EC-Earth3P-VHR, features a triangular truncation at wave number 1279 (hence known as T1279) in spectral space, with a linear N640 reduced Gaussian grid. This corresponds to a spacing of ∼16 km. However, because of the complexity of numerical solutions and parameterizations, the effective resolution (this is the smallest scale that IFS T1279 can fully resolve) is ∼120 km (Abdalla et al., 2013). Vertically, the model features 91 levels, resolving the middle atmosphere up to 0.01 hPa. The model time step during the simulation was 360 s. IFS integrates the revised land surface hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (H-Tessel) model (Balsamo et al., 2009; Hazeleger et al., 2012).

The ocean model is the Nucleus for European Modelling of the Ocean in its version 3.6 (NEMO3.6; Madec, 2008; Madec and the NEMO team, 2016). This is a hydrostatic, finite-difference, free-surface, primitive equation general circulation model. EC-Earth3P-VHR uses the ORCA12 tripolar grid, with the horizontal resolution increasing from the Equator to the poles: ∼9 km at the Equator, ∼7 km at midlatitudes, and ∼2 km near the poles. This corresponds to an effective resolution of ∼45 km (roughly 5 times the ORCA grid spacing; Soufflet et al., 2016). The model uses a z*-coordinate system for the vertical grid and has 75 vertical levels, with the resolution decreasing from 1 m at the surface to 200 m in the deep ocean. The bottom topography is derived from the combination of ETOPO1 (Amante and Eakins, 2009) and GEBCO_08 (Becker et al., 2009). VHR does not include an ocean current feedback (Renault et al., 2023). The sea ice model is the Louvain-la-Neuve sea Ice Model in its version 3 (LIM3) (Vancoppenolle et al., 2012). This is a dynamic–thermodynamic sea ice model, with five ice thickness categories. The time steps are 240 s for NEMO3.6, and 720 s for LIM3 in the EC-Earth3P-VHR.

The atmosphere–land and ocean–sea ice components are coupled through the OASIS (Ocean, Atmosphere, Sea Ice, Soil) coupler, version 3 (OASIS-MCT 3.0) (Valcke and Morel, 2006; Craig et al., 2017). OASIS remaps the atmosphere fluxes onto the ocean grid via nearest-neighbor distance-based Gauss-weighted interpolation. The exchange includes the transfer of momentum, energy, and mass fluxes from the atmosphere to the ocean, while sea surface temperature and sea ice and snow variables are transferred from the ocean to the atmosphere. The remapping of runoff from the atmospheric grid points to runoff areas on the ocean grid re-implemented to be independent of the grid resolution. This was done by introducing an auxiliary model component and relying on the interpolation routines provided by the OASIS coupler. More details on the coupling are provided by Döscher et al. (2022).

EC-Earth3P-VHR (hereafter, VHR) is compared with two lower-resolution global model versions, also run within the PRIMAVERA/HighResMIP project: EC-Earth3P (hereafter, LR; EC-Earth Consortium, 2019) and EC-Earth3P-HR (hereafter, HR; EC-Earth Consortium, 2018). In the atmosphere, they use the T255 (∼107 km) and T511 (∼54.2 km) spectral resolution of the IFS model, respectively (equivalent to an effective resolution of ∼600 and ∼280 km, respectively; Abdalla et al., 2013), both with 91 vertical levels. In the ocean, LR and HR use the ORCA1 (∼100 km) and ORCA025 (∼25 km) tripolar grid, respectively (equivalent to an effective resolution of ∼500 and ∼125 km, respectively; Soufflet et al., 2016), both with 75 vertical levels. They both use the LIM3 sea ice model and the OASIS coupler as well. The LR and HR time steps are respectively 2700 and 900 s in the entire atmosphere, ocean, and sea ice. More details of these two other model versions can be found in Haarsma et al. (2020).

Following the CMIP6 HighResMIP protocol, no additional tuning is applied across resolutions but for a short list of parameters that explicitly change with resolution, particularly for oceanic diffusion and viscosity. The higher resolution in the atmosphere results in a better representation of features such as tropical storms, land–sea transitions, heavy rainfall, and fronts (see Fig. 1 as an example), while in the ocean the increase in resolution allows mesoscale processes to be resolved at a much larger range of latitudes and the representation of finer-resolution bathymetric features and coastlines.

Figure 1 Snapshot of an extratropical storm over the North Atlantic in the winter 1999–2000 in (a) ERA5 and in the (b) VHR, (c) HR, and (d) LR models on their original grids. Shown are daily precipitation rate (mm d⁻¹; blue shading), cloud cover (% of area; gray shading), and sea level pressure (hPa; contours).

2.2 Configuration, workflow setup, and performance optimization

The development and maintenance of the EC-Earth model are supported by the EC-Earth Consortium, which shares model code, configurations, and minimal software infrastructure to operate it. While the LR and HR configurations of EC-Earth-3P were developed in a broad collaboration of all the consortium members participating in PRIMAVERA, VHR's development was primarily completed at the Barcelona Supercomputing Center, in collaboration with the Swedish Meteorological and Hydrological Institute (SMHI) within the ESiWACE2 H2020 project (https://www.esiwace.eu/the-project/past-phases/esiwace2-1/esiwace2, last access: 20 June 2024). The development was conducted on two different supercomputing machines: MareNostrum3 (https://www.bsc.es/marenostrum/marenostrum/mn3, last access: 20 June 2024) and MareNostrum4 (https://www.bsc.es/marenostrum/marenostrum, last access: 20 June 2024). VHR's configuration, at the time of the project, represented one of the most cutting-edge versions of a climate model to run over long timescales. Obtaining a production version of the model, however, entailed (i) generating new grid files; (ii) deploying the initial data; (iii) generating the coupling weights (see below); (iv) creating a new namelist for NEMO; (v) modifying the run scripts to handle the new files and new configuration; (vi) bringing changes from modern versions of the model workflow (Auto-EC-Earth), which, for example, automatizes the call of ELPiN (Tintó et al., 2017; Haarsma et al., 2020) and lets the user fine-tune the distribution of the computational resources in parallel systems; (vii) updating the XIOS (the library for input/output management; https://forge.ipsl.jussieu.fr/ioserver, last access: 30 October 2024) to deal with the land suppression; and (viii) exploring and modifying the configuration parameters to improve the computational throughput of the model execution without losing result accuracy (see below). This presented a significant challenge for both the operations department and the workflow developers, which were required to fine-tune the system to achieve stable runs and minimize the loss of computing hours. Moreover, generating the interpolation weight files to couple the new model grids for the OASIS coupler was particularly challenging. This process could not readily be parallelized at that time in VHR's OASIS3-MCT coupler version (in contrast to more recent ones), and it required collaborating with the OASIS development group. For the workflow, a significant proportion of the effort was devoted to exploiting the hybrid architecture and integrating the dedicated data transfer nodes available in the MareNostrum4 cluster into the workflow software. Additionally, the automatic algorithm that enables the suppression of land grid subdomains in NEMO (ELPiN; Tintó et al., 2017) was incorporated, resulting in a reduction of about 12 % in the required HPC resources (see Haarsma et al., 2020, for more details). Finally, the MareNostrum4 new network (100 Gb Intel Omni-Path Full-Fat Tree), despite its fast and responsive nature, proved to be quite unstable when subjected to high workloads involving multiple concurrent communications, as was the case of the VHR configuration. However, despite the significant challenges, at the end of the ESiWACE2 project (December 2022), the configuration was ready and all the code was versioned and shared with the other partners within the EC-Earth Consortium.

Once deployed, the workflow needed to be made more efficient to be put into operation. Emerging advancements in global climate modeling demand heightened focus on HPC, particularly to accommodate the increasing need for enhanced model resolution (Acosta et al., 2024). An example of such demanding requirements is the VHR configuration, underscoring the need for efficient resource use. In order to address this issue, we conducted a two-fold HPC performance exercise, which involved both a pure computational performance analysis and a scalability study for each model component (IFS and NEMO), complemented with a load balance optimization for the coupling. This analysis concluded that the coupling and output process could be a bottleneck. An optimization was included to package different coupling fields to be sent in the same MPI (Message Passing Interface) communications, reducing the latency and taking advantage of the bandwidth. Additionally, the I/O (input/output) setup was optimized to ensure minimal time was needed to produce the outputs.

While the primary objective of the scalability and load-balance study was to assess the model's efficiency and determine an optimal resource utilization, findings by Acosta et al. (2023) also indicate that enhancing the performance of one component, such as reducing the execution time of IFS, may not necessarily decrease the overall execution time of the coupled model. This discrepancy could stem from a synchronization point at the end of each coupled time step, where both components exchange fields. In cases where other nonoptimized components lag behind, a load rebalance becomes necessary.

Figure 2 Results of the scalability test of the VHR configuration (T1279 IFS and ORCA12 NEMO) at MareNostrum4 (blue line) in simulated years per day (SYPD) for a given number of processors. The orange line shows the ideal case with no loss in computing performance.

Concerning the scalability exercise, we ran a series of tests to balance the resources (computing cores) of the VHR's IFS and NEMO (Fig. 2). To find the most balanced configuration for a given quantity of resources, we followed two different but complementary approaches. The first and most costly one tried to find the optimal distribution by assigning the same number of processors to IFS and NEMO first and moving resources between them alternately; this allowed identifying the intervals for which the model performance increases by using variations of a half-interval search algorithm. The second approach to balance the configuration started from one separate scalability test for each model component that was later used to determine the optimal configuration.

The workflow software Auto-EC-Earth and, by extension, the simulations described here were configured and run with the workflow manager Autosubmit (Manubens-Gil et al., 2016). This Python package facilitates the production of numerical experiments, like the EC-Earth ones, and it allows easily handling experiments with different members, start dates, and initial conditions. The workflow is an oriented graph that includes pre- and post-processing data, the transfer to storage spaces, and the conversion of the output data to CMOR standard, with the details on computing resources needed for each step.

2.3 Simulations

The VHR simulations follow the HighResMIP experimental protocol (Haarsma et al., 2016) and consist of (i) a 50-year spin-up run (spin-up-1950), with initial conditions of temperature and salinity from an ocean state representative of the 1950s (Good et al., 2013, EN4 dataset) and forcing consisting of well-mixed greenhouse gases, including O₃ and aerosol loading for a 1950s (∼10-year mean) climatology; (ii) a 105-year control run (control-1950), starting from the end of spin-up-1950 and keeping the same fixed forcing; (iii) the historical run (hist-1950), starting from the same initial state as the control, but with time-varying external forcing for the period 1950–2014; and (iv) the future scenario run (highres-future), as a continuation of the historical simulation under the CMIP6 SSP5-8.5 scenario (Kriegler et al., 2017) for the period 2015–2050. In this work, VHR's hist-1950 simulation is compared with corresponding hist-1950 runs from LR and HR (Haarsma et al., 2020).

During the model setup, we erroneously applied the EN4 initial conditions at the beginning of all the spin-up runs. While EN4 uses practical salinity and potential temperature, NEMO, which uses the TEOS-10 equation of state, requires absolute salinity and conservative temperature. Nonetheless, the differences between the two temperature and salinity types is indeed small (Pawlowicz, 2013; McDougall et al., 2021), and we expect the error to be minimized throughout the spin-up (see Sect. 3.1).

2.4 Observations and reanalysis

As we mainly aim to evaluate the performance of the EC-Earth3P-VHR configuration and describe the main model biases and characteristics, we focus on the best-observed part of the historical period of the historical simulations, between 1980 and 2014. The three model configurations are compared with the following observational and reanalysis data: near-surface (2 m) air temperature (SAT), zonal winds, sea level pressure, and turbulent fluxes from the ERA5 reanalysis (Hersbach et al., 2020); precipitation rate from the version-2 GPCP dataset (Adler et al., 2003); cloud cover from the version-3 ESA Cloud_cci dataset (ESA CCI-CLOUD; Stengel et al., 2020); potential temperature and salinity of the ocean from the Hadley Centre EN4 (version 4.2.2; Good et al., 2013); sea ice concentration from OSI SAF (OSI-409/OSI-409-a; EUMETSAT Ocean and Sea Ice Satellite Application Facility, 2015); and sea ice volume from GIOMAS (Global Ice–Ocean Modeling and Assimilation System; Zhang and Rothrock, 2003). The period of comparison maximizes data availability and is therefore 1980–2014 for all the cases, except for the GPCP dataset (1983–2014) and the ESA CCI-CLOUD dataset (1982–2014). Biases in sea surface temperature (SST) are very similar to those in SAT and are therefore not shown.

Figure 3 Mean oceanic temperature (in °C) in the LR (yellow), HR (red), and VHR (blue) models in the spin-up runs (0–50-year period), control runs (50–150-year period; solid lines), and historical runs (50–114-year period; dashed lines) in (a) the upper 100 m and (b) the whole ocean. The dashed vertical line marks the end of the spin-up period.

3 Results

3.1 Spin-up phase

Across all three model resolutions, the length of the spin-up (50 years) appears to be insufficient to equilibrate the full ocean (Fig. 3b); in fact, the ocean temperature is still drifting about 0.001–0.002 °C yr⁻¹ (computed over the last 50 years) towards warmer conditions at the end of the control simulation in the three configurations. In the upper ocean, however, VHR shows the smallest warming drift of the three configurations: about 0.00005 °C yr⁻¹ compared to 0.0025 and 0.0062 °C yr⁻¹ in HR and LR, respectively (computed over the last 50 years; Fig. 3a). It is therefore safe to say that an analysis focused on the upper ocean and on the air–sea interface will feature a relatively stable climate in the control simulations. In the historical simulations, the warming of the ocean accelerates due to the CO₂ forcing; after 64 years (year 114 in Fig. 3), the whole ocean warming reaches similar values to those at the end of the control simulations after 100 years in the three model resolutions. Near the surface, the warming trend is much larger. Of the three configurations, VHR is the one with the smallest drift in the control run and the smallest ocean warming in the historical period. Although the three runs start from similar initial conditions derived from an EN4 climatology (Sect. 2.3), VHR is ∼0.4 °C warmer near the surface than LR and HR, especially over the spin-up period. This is likely related to the development of a widespread warm bias over the Southern Ocean (Fig. 4), which we discuss in detail in Sect. 3.6. The trends in global salinity at the end of the control simulations are all smaller than 0.00005 psu yr⁻¹ (computed over the last 50 years; not shown); the three configurations are thus still drifting slightly. As found for the temperature, VHR also shows the smallest drifts of the three configurations (not shown).

Figure 4 Bias in SAT (in K) with respect to ERA5 in the (a) VHR, (b) HR, and (c) LR models for the period 1980–2014. Stippling masks anomalies that are not significant at the 5 % level.

In the following sections, we describe the main characteristics of the VHR compared to LR and HR by focusing on particular regions and biases. This approach should help us highlight the benefits, or lack thereof, due to increased resolution. The main biases in the three model configurations are compared with the observational dataset listed in Sect. 2.4.

3.2 Tropics

A warm bias of 1–2 K is present over the subtropical upwelling regions along the South American and African coasts in the three configurations and shows small variations across them (Fig. 4). The increase in resolution in VHR thus has no clear benefit in reducing it. Past studies have related this bias to an underestimation of the stratocumulus cloud deck (Richter, 2015). This also seems to be the case in the three models, which all show negative cloud biases by about 20 % over all the subtropical upwelling areas, especially along the subtropical Pacific and Atlantic western coasts (Fig. 5). A better-resolved orography near the region does not contribute to reducing the bias either, as suggested in previous studies (Milinski et al., 2016): for example, although VHR shows reduced temperature biases along the Andes compared to HR and LR, it has no effect on the biases over the eastern subtropical Pacific upwelling.

Figure 5 Bias in cloud cover (in %) with respect to ESA CCI-CLOUD (contours in all the panels; in %) in the (a) VHR, (b) HR, and (c) LR models for the period 1982–2014. Stippling masks anomalies that are not significant at the 5 % level.

Figure 6 Bias in precipitation rate (in mm d⁻¹) with respect to GPCP (contours in all the panels; in mm d⁻¹) in the (a) VHR, (b) HR, and (c) LR models for the period 1983–2014. Stippling masks anomalies that are not significant at the 5 % level.

Overall, VHR shows reduced tropical precipitation biases compared to HR and LR (Fig. 6). This is the case, for example, for the double ITCZ bias: this bias is usually characterized by a precipitation excess over the central tropical North Pacific and the western tropical South Pacific and a precipitation deficit over the equatorial Pacific, as LR clearly shows. The dry area over the Equator is reduced with resolution, and the anomaly is even nonsignificant in VHR. This is a clear improvement from increased resolution, and it can be related to a reduced cold bias over the Equator (Fig. 4). In contrast, the precipitation excess over the tropical North Pacific and the Maritime Continent persists into VHR, with only minor reductions of 1–2 mm d⁻¹ compared to HR and LR (Fig. 6). The precipitation excess over the tropical North Pacific suggests a seasonal cycle reaching too far north, while the excess over the Maritime Continent, together with that over the western tropical Atlantic and Indian oceans, suggests an excess in convective precipitation over very warm waters.

Over the tropical Atlantic, the precipitation bias pattern points to an ITCZ anchored to the southwestern part and not reaching the Sahel area. This bias is somewhat reduced in VHR compared to HR and LR, although it is not entirely removed. Over land, the dry bias over northern Brazil, which has been linked to a misrepresentation of the seasonal cycle and extreme events in CMIP6 models (Monteverde et al., 2022), as well as the wet bias along the Andes are not reduced with resolution, either. These positive and negative precipitation biases appear together with positive and negative biases in cloud cover, respectively, related to an overestimation or underestimation in convective clouds (Fig. 5).

3.3 Northern Hemisphere midlatitudes and high latitudes

The largest improvement in the simulated climate from LR to VHR is over the North Atlantic. From south to north, the Gulf Stream representation is much improved in VHR compared to HR and LR, with sharper gradients in temperature and sea surface height (not shown). The position of the Gulf Stream separation is also improved, which leads to a reduction of the warm bias along the US East Coast from LR to VHR (Fig. 4). A paper on a dedicated analysis of the biases over the North Atlantic along the Gulf Stream is currently in preparation (Frigola et al., 2025).

Figure 7 Bias in sea ice thickness (in m; shading) in VHR with respect to GIOMAS for the period 1980–2014. Colored contours are the 15 % value of the sea ice concentration in the LR (orange), HR (red), and VHR (blue) models, as well as in OSI SAF (black) for the period 1980–2014. Panels (a) and (b) are for the Arctic, while panels (c) and (d) are for Antarctica in March (a, c) and September (b, d).

Figure 8 Monthly climatology in the sea ice extent (in 10⁶ km²; a, b) and volume (in 10³ km³; c, d) in the Arctic (a, c) and Antarctica (b, d) in the LR (yellow), HR (red), and VHR (blue) models, as well as in OSI SAF, for sea ice extent and GIOMAS for the volume for the period 1980–2014.

Figure 9 Bias in winter storm track, computed as the standard deviation of the 2–6 d bandpass-filtered daily sea level pressure (in Pa) with respect to ERA5 (contours in all the panels; in Pa) in the (a) VHR, (b) HR, and (c) LR models for the period 1980–2014. Each panel shows anomalies in the boreal winter (DJF; top) and austral winter (JJA; bottom).

Figure 10 Bias in winter zonal wind at 250 hPa (in m s⁻¹) with respect to ERA5 (contours in all the panels; in m s⁻¹) in the (a) VHR, (b) HR, and (c) LR models for the period 1980–2014. Stippling masks anomalies that are not significant at the 5 % level. Each panel shows anomalies in the boreal winter (DJF; top) and austral winter (JJA; bottom).

Farther north, the widespread cold bias up to about 6 K in LR is strongly reduced in HR and even further in VHR, which is the configuration closest to observations (Fig. 4). The cold bias in LR is related to an unrealistically large sea ice extent, which covers the entire Labrador Sea and the western part of the subpolar North Atlantic (Figs. 7 and 8). The reduction of the cold bias between LR and VHR has a deep impact on the climate of the North Atlantic. In the atmosphere aloft, it improves the representation of the boreal winter (DJF) storm track (Fig. 9) and jet (Fig. 10). The boreal winter storm track is overestimated over the subpolar North Atlantic, particularly over the eastern part, in LR, likely related to an excessively strong meridional temperature gradient; by contrast, the VHR storm track is much closer to ERA5 over the North Atlantic.

In the ocean, excessive sea ice leads to a negative salinity bias above 2 psu in the subpolar North Atlantic in LR, which is much reduced in VHR (Fig. 11). Two mechanisms can explain this fresh bias in LR: on the one hand is reduced oceanic salinity transport from subtropical latitudes by a weakened subpolar gyre (not shown) and on the other are errors in the seasonal cycle of the sea ice, during which ice melting would cause an anomalous freshwater input in regions where it is not observed. The negative bias in surface salinity propagates into deeper levels, especially between 300 and 1000 m in the Arctic (Fig. 12). Similarly, the warm subsurface bias at around 40–50° N might also be related to the sea ice excess in the subpolar North Atlantic in LR (Fig. 12). Expanded sea ice in LR causes weaker subpolar gyre strength and associated northward heat transport (not shown), leading to heat accumulation in the intergyre region. However, although this bias is reduced at higher resolutions in HR and VHR, it is still present, suggesting other deficiencies in the formation of intermediate waters in the North Atlantic. The overly large sea ice cover also hampers oceanic deep mixing in the Labrador Sea in LR, whose main regions of deepwater formation are in the Nordic Seas instead (Fig. 13). Oceanic deep mixing takes larger values above 1000 m in VHR and HR in the Labrador Sea. Martin-Martinez et al. (2024) provide a detailed analysis of the characteristics and driving mechanisms of the deepwater formation in the Labrador Sea across the three resolutions and compared to observations.

Figure 11 Sea surface salinity bias (in psu) with respect to EN4 (contours in all the panels; in psu) in the (a) VHR, (b) HR, and (c) LR models for the period 1980–2014. Stippling masks anomalies that are not significant at the 5 % level.

Figure 12 Bias in ocean potential temperature (in K; top) and in salinity (in psu; bottom) with respect to EN4 (contours in all the panels; in K, top, and psu, bottom) in the (a) VHR, (b) HR, and (c) LR models for the period 1980–2014. Stippling masks anomalies that are not significant at the 5 % level. Each panel is separated into the upper and lower 500 m.

Figure 13 Mixed layer depth (in m) in the (a) VHR, (b) HR, and (c) LR models for the period 1980–2014. Northern Hemisphere and Southern Hemisphere values are for March and September, respectively.

Figure 14 Atlantic overturning streamfunction (in Sv) in the (a) VHR, (b) HR, and (c) LR models for the period 1980–2014.

Weak deep mixing results in a relatively weak Atlantic Meridional Overturning Circulation (AMOC; Fig. 14) in LR. The AMOC strength increases with resolution, related to the reduction of the cold bias and sea ice extent bias over the subpolar North Atlantic. The strength of the AMOC in VHR is thus the closest to the observed RAPID strength at 26° N (17±3 Sv, corresponding to the mean and standard deviation, respectively; Frajka-Williams et al., 2019) among the three models: 14±3 Sv in VHR, 12±4 Sv in HR, and 11±2 Sv in LR (computed from monthly streamfunction at 26° N for the period 2004–2014). The structure of the AMOC cell is similar in the three model configurations, with a main positive cell in the upper 3000 m up to 60° N with a maximum at around 30° N, as well as a negative deeper one below with a strength of 2–4 Sv.

In HR, and even more in VHR, the cold bias over the Labrador Sea is replaced by a warm bias (Fig. 4), up to 3–4 K in VHR. This bias also appears in other eddy-rich climate models, related to a stronger ocean heat transport than at lower resolutions in the Atlantic (Roberts et al., 2020b). Over the Nordic Seas, by contrast, a cold bias is present in the three models, although it is somewhat reduced in VHR by 1–2 K compared to LR and HR (Fig. 4). In the three cases, this bias is related to an excessively large sea ice cover in the region (Fig. 7). The warm bias over the Labrador Sea and cold bias over the Nordic Seas in VHR might suggest a misrepresentation of the distribution of oceanic heat transport between the two basins, favoring the westward transport over the northward across-ridge heat transport. It might also or instead be related to a misrepresentation of the sea ice drift across the Denmark Strait (Gutjahr et al., 2022). Relatively weak transport across the strait would lead to an ice deficit in the Labrador Sea (and hence warming) and to ice accumulation in the Nordic Seas (and hence cooling).

On a hemispheric scale, the three models simulate a slightly low Northern Hemisphere sea ice extent, mainly due to the underestimation of the sea ice cover in the Sea of Okhotsk, Baltic Sea, and Labrador Sea in HR and VHR (Fig. 8). By contrast, the three models show an overly large sea ice volume by about 10⁴ km³ compared to GIOMAS (Fig. 9), as they all simulate very thick sea ice in the central Arctic (Fig. 7 for VHR). Anomalously thick ice in the central Arctic would lead to an excess of brine rejection (not shown), which can explain the positive salinity bias above 2 psu in the upper 100–200 m of the Arctic Ocean (Figs. 11 and 12). In VHR, the associated increase in upper-ocean density leads to deeper oceanic mixing than in LR or HR, with a mixed layer depth in the central Arctic that can reach up to 1000 m (Fig. 13).

Over the Pacific, biases tend to be weaker than over the Atlantic. A warm bias of about 1 K develops over the subpolar North Pacific from LR to VHR (Fig. 4), which could explain the negative bias in the boreal winter (DJF) storm track aloft (Fig. 9) and the weaker jet stream over the central Pacific in VHR (Fig. 10).

Over land, the cold bias over the Sahara is reduced with increased resolution (Fig. 4). Similarly, the cold biases over large mountain ranges, such as the Rockies, the Andes, and the Himalaya, up to several degrees in LR are much reduced in VHR (Fig. 4), related to better-resolved orography.

3.4 Southern Ocean

The Southern Ocean is the region where VHR performs the worst compared to HR and LR. The warm bias over the Southern Ocean increases with resolution, up to 4–5 K in VHR compared to 1–2 and 2–3 K for HR and LR, respectively (Fig. 4). It tends to be largest over the Atlantic and Indian sectors of the Southern Ocean and close to the Antarctic coast. Although the warm bias remains generally confined to the upper 100–200 m at around 60° S, it might also be connected to the warm bias at depth between 2000 and 4000 m (Fig. 12).

Two main mechanisms could explain the Southern Ocean warm bias: VHR has the largest cloud cover underestimation of the three models, especially over the Atlantic and Indian sectors, up to 15 % in VHR compared to 5 %–10 % in LR and HR (Fig. 5). Previous studies have related the Southern Ocean warm biases to misrepresentation and underestimation of the mixed-phase clouds, which lead to an excess of shortwave radiation reaching the surface, thereby warming it (e.g., Hwang and Frierson, 2013; Hyder et al., 2018). Connected to the warm bias, VHR also shows the lowest sea ice extent of the three resolutions all year round (Figs. 7 and 8). Although the three models underestimate the Antarctic sea ice extent, in VHR this is nearly half that in observations for the same period (OSI SAF, 1980–2014; EUMETSAT Ocean and Sea Ice Satellite Application Facility, 2015). In terms of sea ice volume (Fig. 8), however, LR shows larger values by about 2×10³ km³ than GIOMAS between November and April, pointing to overly thick sea ice. As for the extent, VHR also shows the lowest sea ice volume, nearly half the values in GIOMAS. The three models show the maximum volume 1 month later than in GIOMAS, in October rather than in September. This contrasts with the Arctic, where the three models capture the general shape of the seasonal cycle.

The surface warming over the Southern Ocean leads to a widespread underestimation of the storm tracks (Fig. 9) and jet stream (Fig. 10) in the austral winter (JJA) in HR and, especially, in VHR compared to LR, which is much closer to ERA5. Although precipitation is also underestimated over the Southern Ocean, especially in VHR, this is not a particularly strong bias, at least compared to those over the tropical regions (Fig. 6).

Late austral summer (September) deep mixing tends to increase by about 200 m from LR to HR and VHR, especially in the Pacific sector. These two latter resolutions show a similar deep mixing mean state, with variations only due to resolution and the better mesoscale representation in VHR (Fig. 13). The underestimation of the storm track over the Southern Ocean therefore does not seem to have an impact on the oceanic mixing below in VHR.

Figure 15 Cross-correlation coefficients between monthly SST and net surface energy flux for the period 1980–2014 in (a) ERA5, (b) J-OFURO3, and the (c) VHR, (d) HR, and (e) LR models. The seasonal cycle and linear trends are removed from the monthly SSTs and energy fluxes before the correlation coefficients are computed. This is done on the original grid in all cases.

3.5 Air–sea coupling

We compare the change in the intensity of air–sea coupling from LR to VHR via the computation of cross-correlation coefficients of the deseasonalized monthly SST and net surface energy flux (Fig. 15). This analysis has extensively been used to study regions in which the ocean tends to drive atmospheric variability (correlation coefficient values approaching 1) or vice versa (correlation coefficient values close to zero; e.g., Bishop et al., 2017; Small et al., 2019). The three model configurations are compared with the ERA5 reanalysis, as done in the previous sections for the biases. To complement the analysis with a non-model-based product, we also include satellite observations of radiative fluxes from J-OFURO3 (Tomita et al., 2019). The two products show overall good agreement, with areas of large correlation coefficient values at the Equator, along the western boundary currents, and over the Southern Ocean (Fig. 15a, b). These areas, nonetheless, tend to be broader in J-OFURO3 than in ERA5.

Over the tropics, the three configurations tend to underestimate the coupling around the Equator, although they all reproduce the band of correlation coefficients of high values along the equatorial Pacific and Atlantic well. However, this band is narrower in LR and HR over the subtropics than it is in ERA5 and J-OFURO3. VHR is thus the closest configuration to the two reference observational products in the region. This result highlights the need for a model resolution finer than 25 km in both the ocean and atmosphere to represent realistic tropical climate interactions, in agreement with the conclusions in Sect. 3.2.

At midlatitudes, the coupling is greatly improved in HR and VHR compared to LR, particularly over the subpolar regions compared to ERA5 and J-OFURO3. LR shows a rather smooth pattern, with very low values in key regions over the Gulf Stream, Kuroshio Current, and Southern Ocean, which suggests that a standard 1° resolution is insufficient to represent a realistic air–sea coupling. VHR and HR show, by contrast, sharper gradients in the correlation coefficient values close to 1 over those regions. This result is consistent with previous studies, which also found a degradation of the air–sea coupling in coarse grids, especially above 1° (e.g., Small et al., 2019). However, VHR shows unrealistic broader areas of higher correlation coefficient values than ERA5 and J-OFURO3 at midlatitudes, degrading the results from HR. One hypothesis for this discrepancy might result from the difference of IFS grid resolution between VHR (T1279) and ERA5 (T639), since the relationship between SST and turbulent fluxes shows certain scale dependency (e.g., Small et al., 2019; Sun and Wu, 2022). However, results do not improve even when regridding VHR onto ERA5 grid before computing the correlation coefficients (not shown). A second hypothesis is the lack of the ocean current feedback in VHR, hence the lack of eddy-killing, which can control the simulated Gulf Stream's dynamics and energy pathways (Renault et al., 2023). However, the pattern of correlation coefficient values remains relatively unchanged when it is computed with a VHR configuration that includes a parameterization that considers the wind adjustment to the ocean current feedback (not shown) (Renault et al., 2019). The results suggest that the VHR's ocean exerts a stronger and more widespread influence on the atmosphere variability than in HR and LR.

Further north, air–sea coupling is overestimated in all the models over the Nordic Seas, likely related to the excess in sea ice in the region and its changes over the seasonal cycle. Together, the results suggest that a realistic air–sea coupling requires grids finer than $\mathrm{1}/\mathrm{4}$° at least, with potential local improvements on a $\mathrm{1}/\mathrm{12}$° grid, especially over the tropics.

4 Discussion and conclusion

This paper presents the eddy-rich configuration of the EC-Earth3P-VHR global model for HighResMIP. We describe both the necessary technical developments to run the model efficiently and the main features of the simulated climate compared to recent observations (1980–2014 period) and to two lower-resolution model configurations (the eddy-present, ∼25 km grid EC-Earth3P-HR; and the non-eddy, ∼100 km grid EC-Earth3P-LR). The EC-Earth3P-VHR (or VHR) uses a comparable atmospheric and oceanic resolution of 10–15 km in a global fully coupled setup, which is, to our knowledge, one of the finest combined grids ever used to date to perform long climate integrations for CMIP (e.g., Small et al., 2014; Chang et al., 2020). Our focus here is on the HighResMIP historical simulation (HighResMIP's hist-1950). This run is part of a larger set of runs, which includes a spin-up run and control runs (HighResMIP's control-1950), a future extension under the SSP8.5 scenario (HighResMIP's highres-future), three hosing simulations forced by idealized Greenland melting, and AMIP sensitivity simulations, all performed within the European PRIMAVERA project and the Spanish STREAM project.

The comparison across the three resolutions (this is VHR, HR, and LR), all with the same physics and no additional tuning, allows identifying regions where increased resolution improves the model performance with respect to observations. One of those regions is the tropics, especially the equatorial Pacific, where the cold tongue bias and the dry bias above are both reduced in VHR compared to HR and LR. Wengel et al. (2021) also report a similar bias reduction in an eddy-resolving configuration of CESM (0.25° resolution in the atmosphere, 0.1° resolution in the ocean), which they link to better-represented mesoscale features, such as tropical instability waves. Similarly, the HadGEM3-GC3.1 global model shows a reduced dry bias over the equatorial Pacific in its configuration with a $\mathrm{1}/\mathrm{12}$° ocean and a 50 km atmosphere (Roberts et al., 2019). By contrast, the eddy-rich MPI-ESM1.2-ER global model ($\mathrm{1}/\mathrm{12}$° ocean as well) shows no evident changes in equatorial precipitation when coupled to a 100 km atmosphere (Gutjahr et al., 2019). Combined, these results suggest that resolutions finer than 25–50 km might be needed in both the atmosphere and ocean to improve surface coupling and reduce biases. However, minimizing equatorial precipitation biases might actually be much more complex than simply increasing model resolution, as found for the ICON global atmosphere–ocean model with a uniform grid spacing of 5 km. Despite its high atmosphere and ocean resolutions, this model still exhibits a strong dry bias over the equatorial Pacific driven by a surface cold bias underneath (Hohenegger et al., 2023; Segura et al., 2022). This model, however, is not directly comparable to other HighResMIP models, as it includes a minimum set of parameterizations. Thus, while convection is directly resolved in ICON, it is parameterized in VHR and the listed models. The incorrect representation of the equatorial SST structure in ICON might instead be related to unresolved subgrid processes (Segura et al., 2022).

The Gulf Stream is another region in which increased model resolution is beneficial, with a reduced temperature biases over the separation region and the central North Atlantic in VHR compared to HR and LR. Such improvements have been related to the resolving of the first baroclinic Rossby radius of deformation over most of the region and/or the exceeding of a critical Reynolds number (e.g., Chassignet and Marshall, 2008) and have been linked to the increase in resolution over the shelf areas to the north of the Gulf Stream (Sein et al., 2017). Similar results have also been reported for the HadGEM3-GC3.1 (Roberts et al., 2019) and MPI-ESM1.2-ER (Gutjahr et al., 2019) global models, both with a $\mathrm{1}/\mathrm{12}$° oceanic grid but coarser atmospheric grids (∼50 and ∼100 km, respectively). This suggests that oceanic resolution is a critical factor for the Gulf Stream representation. Nonetheless, other model features might also be relevant to simulate a realistic Gulf Stream, as no improvement is found in CESM1.3 between a 1° and a 0.1° oceanic grid, for which the Gulf Stream separation occurs too far north (Chang et al., 2020). One of the many potential reasons behind the discrepancy might be the obvious difference in the number of atmospheric vertical levels: 91 in VHR, 85 in HadGEM3-GC3.1 (Roberts et al., 2019), 95 in MPI-ESM1.2-ER (Gutjahr et al., 2019), and only 30 in CESM1.3 (Meehl et al., 2019), which is expected to degrade the representation of key stratosphere–troposphere interactions affecting North Atlantic variability, and, by extension, the wind field, which is critical for the Gulf Stream separation. As nicely summarized in Chassignet and Marshall (2008), however, “The Gulf Stream separation, indeed, turns out to be quite sensitive to a variety of other factors such as subgrid scale parametrization, subpolar gyre strength and water mass properties, [deep western boundary current] strength, representation of topography, and the choice of model grid”. A realistic representation of the Gulf Stream is crucial for the North Atlantic and European climate. SST biases in the Gulf Stream can drive not only local changes over the North Atlantic, but also a large-scale dynamic response over remote regions of the Northern Hemisphere through a quasi-zonal planetary barotropic Rossby wave response (Lee et al., 2018). Similarly, a more realistic, farther-south Gulf Stream has been shown to shift north in simulations with increased CO₂ in models at eddy-rich resolutions (Saba et al., 2016; Moreno-Chamarro et al., 2021). This shift would lead to amplified warming of the US East Coast region, which might be consistent with the anomalous warming observed in the Gulf Stream area in recent decades (Pershing et al., 2015; Todd and Ren, 2023). Reducing biases in the Gulf Stream area is therefore key to reproducing a realistic atmospheric circulation and to the sensitivity of the response to an external forcing.

Mainly related to increased atmospheric resolution, VHR also shows reduced precipitation biases over mountain ranges all over the world. This suggests that VHR might provide more realistic regional information on precipitation variability and future changes than lower-resolution models can. Giorgi et al. (2016), in fact, showed that increased model resolution leads to stronger summer precipitation changes over the Alpine region using climate change projections from a regional atmospheric model with a ∼12 km grid. VHR uses a similar resolution but on a global scale, without the need to be constrained by lower-resolution models.

On the negative side, we find that increased model resolution alone can be insufficient to reduce important and well-known biases in the climate or even cause model degradation in VHR. The warm bias over the coastal tropical upwelling areas, the Southern Ocean warm bias, and the rainfall excess bias over warm tropical waters all persist or even increase in VHR compared to HR and LR. These biases point to deficiencies in the model physics, especially in the atmosphere, and more particularly in the cloud parameterizations. In VHR, the warm bias over both eastern tropical upwelling areas and the Southern Ocean is connected to negative biases in cloud cover. This reinforces the established idea that insufficient stratocumulus decks over the upwelling areas (e.g., Richter, 2015) and mixed-phase clouds over the Southern Ocean (e.g., Hyder et al., 2018) play key roles in setting up those biases. Cloud biases can be particularly insensitive to increases in model resolution, both in the ocean and atmosphere, from ∼100 km grids to 25–50 km grids (Moreno-Chamarro et al., 2022). Yet, for example, improved cloud microphysics closer to observations have been shown to help reduce shortwave radiation biases over the Southern Ocean in the Met Office's Unified Model (Varma et al., 2020). Reducing these biases as much as possible is critical, since they can have wider, global impacts on the climate, driving, for example, additional biases in tropical precipitation through the effect on the global energy budget (e.g., Hwang and Frierson, 2013; Hawcroft et al., 2017).

It is interesting to note, nonetheless, that although LR, HR, and VHR all share the same cloud scheme, it is VHR that develops the strongest Southern Ocean bias. This might be related to the lack of additional model tuning from LR to HR and VHR. Rackow et al. (2024) showed that tuning the top-of-the-atmosphere radiation contributed to reducing the warming excess over the Southern Ocean in the IFS-FESOM global model at ∼5 km resolution. The HighResMIP protocol suggests that no tuning be performed across resolutions to ensure that any changes in the simulated climate can solely be attributed to changes in resolution (Haarsma et al., 2016). This approach can lead to undesired model degradation: for example, the untuned, low-resolution ECMWF model for HighResMIP shows an overly weak AMOC and a large cold bias over the North Atlantic compared to its well-tuned, high-resolution counterpart (Roberts et al., 2018a). This can hinder model comparison and a clear understanding of the effect of model resolution, as biases can have large-scale climatic impacts (e.g., Hwang and Frierson, 2013; Hawcroft et al., 2017; Lee et al., 2018) and affect the response sensitivity to forcing (e.g., McGee et al., 2018).

With respect to the spin-up, the HighResMIP protocol suggests a 50-year period (Haarsma et al., 2016). For all the configurations, this period is insufficient to equilibrate the full ocean, although the upper 100 m equilibrates faster than the lower part, and VHR does it faster and appears more stable after 100 years than HR and LR. The eddy-rich HadGEM3-GC3.1 also shows smaller drifts at the end of the 50-year period than its lowest-resolution versions (Roberts et al., 2019). By contrast, for CESM1.3, the low- and high-resolution configurations only show a more stable climate after 150 years, related to a strong top-of-the-atmosphere energy imbalance (Chang et al., 2020). This led the authors to propose “150 to 200 years of model spin-up as a future strategy for initializing HR climate model simulations” (Chang et al., 2020). However, considering how computationally expensive these simulations are, new techniques might need to be introduced to tune and spin these models up faster and for longer. As much as tuning can still be “artisanal in character” at many research centers (Mauritsen et al., 2012), new and faster methods are being implemented to speed up the exploration of the space of parameters to find the best fit with observations. These methods include, for example, machine learning (Hourdin et al., 2021), simplified configurations (Wan et al., 2014), adjoints (Lyu et al., 2018), or model emulators (Williamson et al., 2013). Additional techniques have also been proposed to spin models up faster with much less computational cost; these include using, for example, Newton–Krylov methods (Bernsen et al., 2010; Merlis and Khatiwala, 2008) or replacing the atmosphere model with model data (Lofverstrom et al., 2020). Implementing similar techniques in future HR and VHR simulations would help accelerate both the spin-up and tuning phases.

To summarize, we present the eddy-rich version of the EC-Earth global climate model, EC-Earth3P-VHR, with atmospheric and oceanic resolutions of 10–15 km. The analysis of its main climate features reveals improvements with respect to two lower-resolution versions, such as a reduced dry equatorial bias over the Pacific, a more realistic Gulf Stream representation, and more accurate rainfall over mountain areas. Other biases persist or degrade, such as the warm biases over the subtropical upwelling regions and Southern Ocean, the tropical precipitation excess, and the excess in sea ice volume and oceanic deep mixing in the Arctic. VHR's global resolution is at a similar level as many regional models, such as those participating in CORDEX, and it is much finer than most of the standard CMIP models. This opens a window of opportunity for model comparison and evaluation, as well as process understanding of a much more realistic present-day and future climate on a more regional scale.