3.1 Scalability

The results of the scalability analyses of the atmosphere (IFS) and ocean (NEMO) components of EC-Earth3P-HR are shown in Fig. 1 and for the fully coupled model in Fig. 2. Acosta et al. (2016) showed that, while for coupled application the load balance between components has to be taken into account in the scalability process, the process needs to start with a scalability analysis of each individual component. Moreover, the user could experience that the speeding up of one component (e.g., the reduction of the execution time of IFS) does not reduce the execution time of the coupled application. This could be because there is one synchronization point at the end each coupled time step where both components exchange fields. If the other non-optimized components are slower, a load rebalance will be required. The final choice depends on the specific problem, where either time or energy can be minimized. In Sect. 3.2, we describe how the optimal load balance between the two components, where NEMO is the slowest component, was achieved (Acosta et al., 2016).

3.2 Bottlenecks

For the performance analysis, the individual model components (IFS, NEMO and OASIS) are benchmarked and analyzed using a methodology based on extracting traces from real executions. These traces are displayed using the PARAVER (PARAllel Visualization and Events Representation) software and processed to discover possible bottlenecks (Acosta et al., 2016). Eliminating these bottlenecks not only involves an adjustment of the model configuration and a balance of the number of cores devoted to each one of its components but also modifications of the code itself and work on the parallel programming model adopted in the different components.

Figure 3(a) PARAVER view of the NEMO and IFS components in an EC-Earth3P-HR model execution for two time steps including the coupling process. The horizontal lines give the behavior of the different processes (1 to 512 for IFS and 513 to 536 for NEMO) as a function of time. Each color corresponds to a different MPI communication function. See text for explanation. Panel (b) is the same as (a) but when optimization options “opt” and “gathering” for coupling are activated.

Download

The first step of a performance analysis consists in analyzing parallel programming model codes using targeted performance tools. Figure 3a illustrates an example of the performance tool's output from one single EC-Earth3P-HR model execution as provided by the PARAVER tool, focusing only on its two main components: NEMO and IFS. This figure is very useful for determining the communications within the model and identify sources of bottlenecks, especially those resulting from communication between components. It displays the communication pattern as a function of time. The vertical axis corresponds to the different processes executing the model, the top part for IFS and the lower part for NEMO. The different colors correspond to different MPI communication functions, except the light blue, which corresponds to no communication. Red, yellow and purple colors are related to MPI communications. The green color represents the waiting time needed to synchronize the coupled model for the next time step, which means an unloaded balance in the execution. In summary, light blue areas are pure computation and should be maximized. On the other hand, yellow, red and purple are representing overhead from parallel computation and should be minimized if possible. Additionally, green areas are preferably to be also reduced, for example, increasing the number of parallel resources of the slowest component, but no optimizations are needed. From this analysis, several things can be concluded related to the overhead from parallel computation:

1. Figure 3 shows the coupling cost from a computational point of view, including one regular time step of IFS and NEMO and one time step including the coupling process. In the top part of Fig. 3a, we notice that during the first half of the first time step, the IFS component model reserves most of its processors for execution (512 processes). To simplify, it can be said that the first half of the time step has less MPI communication, with more computation-only regions, while the second half of the time step is primarily about broadcasting messages (yellow and white color block), which corresponds to the coupling computation and to send/receive files from the atmospheric to the ocean model. These calculations impact the scalability of the code dramatically. This configuration increases the overhead when more and more processes are used and represents more than 50 % of time execution when 1024 processes are used. The coupling process can be analyzed in detail in Fig. 3a (coupling zoom, top image), where the same pattern of communications is repeated four times. This occurs because the different fields from IFS to NEMO are sent in three different groups, followed by an additional group of fields sent from IFS to the runoff mapper component. The communication of three different groups of fields to the same component is not taking advantage of the bandwidth of the network, thus increasing the overhead produced by MPI communications. However, these three groups are using the same interpolation method and they could be gathered into the same group.

2. From other parts of the application (not shown in the figure), we also notice the expensive cost of the IFS output process for each time step. A master process gathers the data from all MPI subdomains and prints the complete outputs at a regular time interval of 3–6 h. During this process, the rest of processes are waiting for this step to be completed. Due to the large data volumes, this sequential process is very costly, increasing the execution time of IFS by about 30 % when outputs are required, compared to the regular time step of IFS (without output).

3. The bottom part of Fig. 3a shows that the communication in NEMO is not very effective and that a large part of it is devoted to global communications, which appear in purple. Those communications belong to the horizontal diffusion routine, inside the ice model (LIM3) used in NEMO. The high frequency of communications in this routine prevented the model from scaling. More information about MPI overhead of NEMO can be found in Tintó et al. (2019).

4. Due to the domain decomposition used by NEMO, some of the MPI processes, which are used to run part of the ocean domain in parallel, were computing without use. This is because domain decomposition is done on a regular grid and a mask is used to discriminate between land and sea points. The mask creates subdomains of land points whose calculations are not used. This is illustrated in Fig. 4, which shows a particular case in which 12 % of the depicted subdomains do not contain any sea point.

Figure 4Domain decomposition of a tripolar grid of the ORCA family with a resolution of 1.0^∘ into 128 subdomains (16×8). Subdomains marked with a black dot do not contain any ocean grid point.

3.3 New optimizations for the specific configuration

According to the profiling analysis done, different optimizations were implemented to improve the computational efficiency of the model:

1. The optimization (“opt”) option of OASIS3-MCT was used. This activates an optimized global conservation transformation. Using this option, the coupling time from IFS to NEMO is reduced by 90 % for EC-Earth3P-HR. This is because all-to-one/one-to-all MPI communications are replaced by global communications (gather/scatter and reduction) and the coupling calculations are done by all the IFS processes instead of only the IFS master process.

   Another functionality of OASIS consists in gathering all fields sent from IFS to NEMO in a single group (Acosta et al., 2016). Coupling field gathering, an option offered by OASIS3-MCT, can be used to optimize coupling exchanges between components. The results show that gathering all the fields that use similar coupling transformations reduces the coupling overhead. This happens because OASIS3-MCT is able to communicate and interpolate all of the fields gathered at the same time. Figure 3a (coupling zoom, bottom image) proves that the collection of the first three groups reduces the communication patterns from four to two, where the execution time of this part is reduced significantly (40 %).

   Figure 3b shows the execution when “opt” and “gathering” options are used, with the 90 % reduction in coupling time clearly visible (large green section). In the case of the first time step in the trace, the coupling time is replaced by waiting time, since NEMO is finishing its time step and both components have to exchange fields at the end of the time step.

2. For the output problem, the integration of XIOS as the I/O server for all components of EC-Earth can increase performance dramatically. XIOS is already used for the ocean component NEMO and the I/O server receiving also all the data from IFS processes and doing the output work in parallel and in an asynchronous way is the best solution to remove the sequential process when an IFS master process is required to do this work. This is being developed and will be included in the next version of EC-Earth.

3. Based on the performance analysis, the amount of MPI communications can be reduced (Tintó et al., 2019), achieving a significant improvement in the maximum model throughput. In the case of EC-Earth3P-HR, this translated into a reduction of 46 % in the final execution time.

4. Using the tool ELPiN (Exclude Land Processes in NEMO), the optimal domain decomposition for NEMO has been implemented (Tintó et al., 2017), with computation of only ocean subdomains and finding the most efficient number of MPI processes. This substantially improves both the throughput and the efficiency (in the case of 2048 processor cores, 41 % faster using 25 % less resources). The increase in throughput was due to fewer computations and, related to that, fewer communications. In addition, ELPiN allows for the optimal use of the available resources in the domain decomposition depending on the shape and overlap of the subdomains.

3.4 Post-processing and data output

At the T511L91 resolution, the HighResMIP data request translates into an unprecedented data volume for EC-Earth. Because the atmosphere component (IFS) is originally a numerical weather prediction (NWP) model, it contains no built-in functionality for time averaging the data stream during the simulation. The model was therefore configured to produce the requested three-dimensional fields (except radiative fluxes on model levels, which cannot be output by the IFS) on a 6-hourly basis and surface fields with 3-hourly frequency. As a consequence, the final daily and monthly averages for instantaneous fields have been produced from sampling at these frequencies, whereas fluxes are accumulated in the IFS at every time step. Vertical interpolation to requested pressure or height levels is performed by the model itself.

For the ocean model, the post-processing is done within NEMO by the XIOS library which can launch multiple processes writing NetCDF files in parallel, alleviating the I/O footprint during the model run. The XIOS configuration XML files were extended to produce as many of the ocean and sea-ice variables as possible.

The combination of the large raw model output volume, the increased complexity of the requested data and the new format of the CMOR tables (Climate Model Output Rewriter, an output format in conformance with all the CMIP standards) required a major revision of the existing post-processing software. This has resulted in the development of the ece2cmor3 package. It is a Python package that uses Climate Data Operators (CDO) (CDO, 2015) bindings for (i) selecting variables and vertical levels, (ii) time averaging (or taking daily extrema), (iii) mapping the spectral and grid-point atmospheric fields to a regular Gaussian grid and (iv) computing derived variables by some arithmetic combination of the original model fields. Finally, ece2cmor3 uses the PCMDI CMOR library for the production of NetCDF files with the appropriate format and metadata. The latter is the only supported step for the ocean output.

To speed up the atmosphere post-processing, the tool can run multiple CDO commands in parallel for various requested variables. Furthermore, we optimized the ordering of operations, performing the expensive spectral transforms on time-averaged fields wherever possible. We also point out that the entire procedure is driven by the data request; i.e., all post-processing operations are set up by parsing the CMOR tables and a single dictionary relating EC-Earth variables and CMOR variables. This should make the software easy to maintain with respect to changes in the data request and hence useful for future CMIP6 experiments.

4.1 Outline of HighResMIP protocol

The protocol of the HighResMIP simulations consists of tiers 1, 2 and 3 experiments that represent simulations of different priority (1 being the highest, 3 the lowest) and a spin-up procedure. The protocol also excludes specific tuning for the high-resolution version compared to the standard-resolution version. Below, we give a short summary of the protocol. The experiment names in the CMIP6 database are given in italics.

• Tier 1: forced-atmosphere simulations 1950–2014; highresSST-present. The Tier 1 experiments are atmosphere-only simulations forced using observed sea surface temperature for the period 1950–2014.

• Tier 2: coupled simulations (1950–2050). The period of the coupled simulations is restricted to 100 years because of the computational burden brought about by the model resolution and the limited computer resources. The period 1950–2050 covers historical multi-decadal variability and near-term climate change. The coupled simulations consist of a spin-up, control, historical and future simulations.

  • Spin-up simulation; spinup-1950. Due to the large computer resources needed, a long spin-up to (near) complete equilibrium is not possible at high resolution. Therefore, as an alternative approach, an analyzed ocean state representative of the 1950s is used as the initial condition for temperature and salinity (Good et al., 2013, EN4 data set). To reduce the large initial drift, a spin-up of about 50 years is made using constant 1950s forcing. The forcing consists of greenhouse gases (GHGs), including O₃ and aerosol loading for a 1950s (∼ 10-year mean) climatology. Output from the initial 50-year spin-up is saved to enable analysis of multi-model drift and bias, something that was not possible in previous CMIP exercises, with the potential to better understand the processes causing drift in different models.

  • Control simulation; control-1950. This is the HighResMIP equivalent of the pre-industrial control but uses fixed 1950s forcing. The length of the control simulation should be at least as long as the historical plus future transient simulations. The initial state is obtained from the spin-up simulation.

  • Historical simulation; hist-1950. This is the coupled historical simulation for the period 1950–2014, using the same initial state from the spin-up as the control run.

  • Future simulation; highres-future. This is the coupled scenario simulation 2015–2050, effectively a continuation of the hist-1950 experiment into the future. For the future period, the forcing fields are based on the CMIP6 SSP5-8.5 scenario.

• Tier 3: forced-atmosphere 2015–2050 (2100); highresSST-future. The Tier 3 simulation is an extension of the Tier 1 atmosphere-only simulation to 2050, with an option to continue to 2100. To allow comparison with the coupled integrations, the same scenario forcing as for Tier 2 (SSP5-8.5) is used.

A schematic representation of the HighResMIP simulations is given in Fig. 5.

Figure 5Schematic representation of the HighResMIP simulations.

Download

Table 2Overview of the HighResMIP simulations of EC-Earth3P-HR and EC-Earth3P.

Download Print Version | Download XLSX

4.2 Main results of EC-Earth3P(-HR) HighResMIP simulations

For each of the HighResMIP tiers, more than one simulation was produced. An overview of the simulations is given in Table 2.

The data are stored on the JASMIN server at the Centre for Environmental Data Analysis (CEDA) (https://www.ceda.ac.uk/services/jasmin/, last access: 29 July 2020) and available from Earth System Grid Federation (ESGF). During the PRIMAVERA project, the data were analyzed on the JASMIN server. For the highresSST-present and highresSST-future simulations, the ensemble members were started from perturbed initial states. These were created by adding small random perturbations from a uniform distribution over [$-\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$, $+\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$] degree to the three-dimensional temperature field. For the control-1950 and hist-1950, the end of the spin-up was taken as the initial condition of the first member. For the two extra members, the initial conditions were generated by continuing the spin-up for 5 years after perturbing the fields that are exchanged between atmosphere and ocean. The highres-future members are the continuation of the hist-1950 members.

The Atlantic meridional overturning circulation (AMOC) in the control-1950 of EC-Earth3P had unrealistically low values of less than 10 Sv. It was therefore decided to change the ocean mixing parameters, which improved the AMOC. The main difference compared to the first ensemble member of EC-Earth3P is that the parameterization of the penetration of turbulent kinetic energy (TKE) below the mixed layer due to internal and inertial waves is switched off (nn_etau = 0; Madec and the NEMO team, 2016). The mixing below the mixed layer is an ad hoc parameterization into the TKE scheme (Rodgers et al. 2014) and is meant to account for observed processes that affect the density structure of the ocean's boundary layer. In EC-Earth3P, this penetration of TKE below the mixed layer caused a too-deep surface layer of warm summer water masses in the North Atlantic convection areas, which led to a breakdown of the Labrador Sea convection within a few years and a strongly underestimated AMOC in EC-Earth3P. An additional minor modification compared to ensemble member 1 is an increased tuning parameter rn_lc (= 0.2) in the TKE turbulent closure scheme that directly relates to the vertical velocity profile of the Langmuir cell circulation. Consequently, the Langmuir cell circulation is strengthened.

The new mixing scheme was also applied to EC-Earth3P-HR to ensure the same set of parameters for both versions of EC-Earth3P(-HR). The simulations with the new ocean mixing are denoted with “p2” for the coupled simulations in Table 2. The atmosphere is unchanged, and therefore the atmosphere simulations are denoted as “p1”. Because of the unrealistically low AMOC in EC-Earth3P in the “p1” simulations, we focus on “p2” for the coupled simulations.

Below, we will briefly discuss the mean climate and variability of the highresSST-present, control-1950 and hist-1950 simulations. The main differences between EC-Earth3P and EC-Earth3P-HR will be highlighted. In addition, the spin-up procedure for the coupled simulations, spinup-1950, will be outlined. A more extensive analysis of the HighResMIP simulations will be presented in forthcoming papers.

4.2.1  highresSST-present

The highresSST-present simulations will be compared with ERA5 (Hersbach et al., 2020) (1979–2014), except for precipitation, where Global Precipitation Climatology Project v2.3 (1979–2014) (Adler et al., 2003) data will be used. EC-Earth, GPCP and ERA5 data are regridded to a common grid (N128) before comparison. Seasonal means (December–February (DJF) and June–August (JJA)) will be analyzed. Ensemble mean fields will be displayed.

Figure 6SAT: bias (^∘C) of EC-Earth3P-HR with respect to ERA5 for the period 1979–2014: (a) DJF and (b) JJA. Global means of SAT for EC-Earth3P-HR are 11.01 (DJF) and 15.85 (JJA), and for ERA5 12.43 (DJF) and 15.95 (JJA). RMSEs of EC-Earth3P-HR with respect to ERA5 are 1.25 (DJF) and 1.06 (JJA).

Figure 7MSLP: bias (hPa) of EC-Earth3P-HR with respect to ERA5 for the period 1979–2014: (a) DJF and (b) JJA. Global means of MSLP for EC-Earth3P-HR are 1011.3 (DJF) and 1009.4 (JJA), and for ERA5 1011.53 (DJF) and 1011.24 (JJA). RMSEs of EC-Earth3P-HR with respect to ERA5 are 1.11 (DJF) and 1.27 (JJA).

Figure 8Precipitation and bias of EC-Earth3P-HR with respect to GPCP (mm d⁻¹) for the period 1979–2014: (a, c) DJF and (b, d) JJA. Global means of precipitation for EC-Earth3P-HR are 2.91 (DJF) and 3.25 (JJA), and for ERA5 2.70 (DJF) and 2.71 (JJA). RMSEs of EC-Earth3P-HR with respect to ERA5 are 1.06 (DJF) and 1.44 (JJA).

Due to the prescribed SST, the largest surface air temperature (SAT) biases are over the continents (Fig. 6). The largest negative biases are over central Africa for DJF and Alaska in JJA, while the largest positive biases are located over Antarctica in JJA and northeastern Siberia in DJF. Over most areas, EC-Earth3P-HR is slightly too cold. Over most of the tropics, the mean sea level pressure (MSLP) is underestimated, whereas over Antarctica and the surrounding regions of the Southern Ocean it has a strong positive bias (Fig. 7). Also noteworthy is the positive bias south of Greenland during DJF. The largest precipitation errors are seen in the tropics over the warm pool regions in the Pacific and the Atlantic with too much precipitation (Fig. 8). The planetary wave structure of the geopotential height at 500 hPa (Z500) during DJF is well represented, with the exception of the region south of Greenland (Fig. 9), which is consistent with the MSLP bias (Fig. 7a). The physical causes of the aforementioned biases can include a wide range of deficiencies in the parameterizations of cloud physics, land surface and snow, to mention a few. In forthcoming papers, this will be investigated in further detail.

Figure 9(a) Stationary eddy component (departure from zonal mean) of EC-Earth3P-HR of the 500 hPa geopotential height (m) in boreal winter; (b) the difference with ERA5. Note the difference in color scale between the two panels.

Figure 10Differences between EC-Earth3P-HR and EC-Earth3P for SAT (^∘C) (a, b), MSLP (hPa) (c, d) and precipitation (mm d⁻¹) (e, f) for DJF (a, c, e) and JJA (b, d, f).

Doubling of the atmospheric horizontal resolution has only a modest impact on the large-scale structures of the main meteorological variables, as illustrated by the global MSLP, SAT and precipitation (Fig. 10). For SAT, the differences are generally less than 1 K; for MSLP, they are 1 hPa, except in the polar regions. A remarkable result is the worsening of the bias over Antarctica during JJA. Because the dynamics of the polar vortex, which is sensitive to horizontal resolution, is strongest during austral winter, we speculate that this enhanced bias is associated with it. The exact mechanism falls outside the scope of this basic validation and will be explored in forthcoming studies. For precipitation, the difference can be larger than 1.5 mm d⁻¹ in the tropics. It is possible to conclude that the increase of resolution does not have a clear positive impact on the climatology of any of those variables. For instance, for precipitation, it results in an increase of the wet bias over the warm pool (compare with Fig. 8). Also measured by the root mean square error (RMSE) (see figure captions for the numbers), the impact of resolution is small, on the order of 10 % or less depending on the variable and the season. Enhancing resolution reduces the RMSE for SAT and MSLP, whereas it slightly increases for precipitation (from 1.04 to 1.06 in DJF and from 1.35 to 1.44 in JJA).

4.2.2  spinup-1950

As discussed in the outline of the HighResMIP protocol, the spin-up was started from an initial state that is based on observations for 1950. For the ocean, this is the EN4 ocean reanalysis (Good et al., 2013) averaged over the 1950–1954 period, with 3 m sea-ice thickness in the Arctic and 1 m in the Antarctic. The atmosphere–land system was initialized from ERA-20C for 1 January 1950 and spun up for 20 years to let the soil moisture reach equilibrium. For the ocean, no data assimilation has been performed, which can result in imbalances between the density and velocity fields, giving rise to initial shocks and waves.

During the first years of the spin-up, there is a strong drift in the model climate (not shown). For the fast components of the climate system like the atmosphere and the mixed layer of the ocean, the adjustment is on the order of 1 year, whereas the slow components such as the deep ocean require a thousand years or more to reach equilibrium. For the land component, this is on the order of a decade. As a consequence, after a spin-up of 50 years, the atmosphere, land and upper ocean are approximately in equilibrium, while the deeper ocean is still drifting. The largest drift occurs in the layer of 100–1000 m with a drift of 0.5 ^∘C per century. This drift also has an impact on the fast components of the climate system, which therefore still might reveal trends.

Figure 11(a) Global mean averaged annual SAT (^∘C) in control-1950 for the three members of EC-Earth3P (red colors) and EC-Earth3P-HR (grey colors). (b) Global mean averaged net surface heat flux (W m⁻²) in control-1950 of EC-Earth3P (red) and EC-Earth3P-HR (black), displayed only for one member (r1i1p2f1) of each model for clarity; other members display similar behavior.

Download

4.2.3  control-1950

After the spin-up, the SAT each of the three members of EC-Earth3P-HR is in quasi-equilibrium and the global mean temperature oscillates around 13.9 ^∘C (Fig. 11a, black). The ocean is still warming, as expressed by a negative net surface heat flux on the order of −1.5 W m⁻² (positive is upward) (Fig. 11b, black). This imbalance is reduced during the simulation but without an indication that the model is getting close to its equilibrium state.

Figure 12Ensemble mean SAT (^∘C) of the averaged last 10 years (2040–2049) minus the averaged first 10 years (1950–1959) of the control-1950 simulations of EC-Earth3P.

Contrary to EC-Earth3P-HR, the global annual mean SAT of EC-Earth3P displays a significant upward trend, with an indication of stabilizing after about 35 years (Fig. 11a, red). This warming trend is caused by a large warming of the North Atlantic, as revealed in Fig. 12, showing the difference between the first and last 10 years of the control-1950 run. This warming is caused by the activation of the deep convection in the Labrador Sea (not shown) that started about 10 years after the beginning of the control simulation, which was absent in the spin-up run. Associated with that, the AMOC also shows an upward trend (see Fig. 17 below). This switch to a warmer state does not strongly affect the slow warming of the deeper ocean, which is reflected in a similar behavior of the net surface heat flux to that for EC-Earth3P-HR (Fig. 11b). The reasons for the initial absence of deep convection in the Labrador Sea in EC-Earth3P and the difference with EC-Earth3P-HR are not clear and presently under investigation. Possible candidates are that the differences in ocean resolution affect the sea-ice dynamics and deep convection, but also changes in ocean temperature and salinity distribution may play a role.

The control-1950 experiment is also analyzed to evaluate model performance of internally generated variability in the coupled system; the targets are the El Niño–Southern Oscillation (ENSO), the North Atlantic Oscillation (NAO), sudden stratospheric warmings (SSWs) and the Atlantic meridional overturning circulation (AMOC).

Figure 13Monthly standard deviation of the NINO3.4 SST index: EC-Earth3P (red) EC-Earth3P-HR (blue) from control-1950 and detrended Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST) over 1900–2010 (black).

Download

Figure 14(a–c) Boreal winter SST standard deviation from control-1950 in EC-Earth3P (a), EC-Earth3P-HR (b) and detrended HadISST (c); overplotted with contours are the corresponding climatology (contour interval 2 ^∘C). Bottom: regression of SST anomalies onto the NINO3.4 index from control-1950 in EC-Earth3P (d), EC-Earth3P-HR (e) and detrended HadISST (f); overplotted with contours are the corresponding regression of 500 hPa geopotential height anomalies (c.i. 2.5 m), with ERA-Interim in panel (f). The observational period is 1979–2014. (g) Difference in SST standard deviation between EC-Earth3P-HR (b) and EC-Earth3P (a).

SSWs

Also, the simulation of SSW occurrence is assessed (Fig. 16); the identification follows the criterion in Palmeiro et al. (2015). The decadal frequency of SSWs in EC-Earth is about eight events per decade regardless of model resolution, which is underestimated when compared to ERA-Interim (∼ 11 events per decade) but in the range of observational uncertainty (e.g., Palmeiro et al., 2015; Ayarzagüena et al., 2019). The same underestimation was diagnosed in EC-Earth3.1 (Palmeiro et al., 2020a). The reduced amount of SSWs is probably associated with a too-strong bias at the core of the polar vortex, still present in EC-Earth3.3 (Palmeiro et al., 2020b). It is thus concluded that increasing horizontal resolution does not affect the model bias in the strength of the polar vortex. The seasonal cycle of SSWs in reanalysis is quite robust over the satellite period, showing one maximum in December–January and another one in February–March (Ayarzagüena et al., 2019), which was properly captured by EC-Earth3.1 in the control, coupled simulations with fixed radiative forcing at the year 2000 (Palmeiro et al., 2020a). Here, in control-1950, EC-Earth does not reproduce such a bimodal cycle, with EC-Earth3P-HR (blue) yielding a peak in January–February and EC-Earth3P (red) two relative maxima in January and March. Interestingly, the seasonal cycle of SSWs over the historical, pre-satellite period shows a different distribution with a prominent maximum in mid-winter and a secondary peak in late winter, although it is less robust among reanalysis products (Ayarzagüena et al., 2019). The impact of the radiative forcing on SSW occurrence deserves further research.

Figure 15(a–c) Regression of 500 hPa geopotential height anomalies from control-1950 in EC-Earth3P (a), EC-Earth3P-HR (b) and detrended ERA-Interim (c) onto the corresponding leading principal component, i.e., the NAO index. (d–f) Leading EOF of winter SLP anomalies over the North Atlantic–European region 20–90^∘ N, 90^∘ W–40^∘ E, from control-1950 in EC-Earth3P (d), EC-Earth3P-HR (e) and detrended ERA-Interim (f); the corresponding fraction of explained variance is indicated in the title.

4.2.4  hist-1950

The hist-1950 ensemble simulations differ from the control-1950 simulations in terms of the historical GHG and aerosol concentrations. The global mean annual temperature in EC-Earth3P-HR displays an increase similar to the ERA5 data set (Fig. 18a). The warming seems to be slightly larger in the model. We remind the reader, however, of the enhanced observed warming after 2014, which might result in a similar trend in the model simulations compared to observations up to present day. The cooling due to the Mt. Pinatubo eruption in 1991 is clearly visible in all members and the ensemble mean. The amplitude and period compare well with ERA5. On its part, the AMOC in EC-Earth3P-HR reveals a clear downward trend in particular from the 1990s onward (Fig. 17b, black). This is consistent with a slowdown of the Atlantic overturning due to global warming in CMIP5 models (Cheng et al., 2013).

Similarly to control-1950, the hist-1950 simulations with EC-Earth3P show an upward drift in SAT (Fig. 18a, red) and AMOC (Fig. 17b, red) that are smaller (SAT) or absent (AMOC) in EC-Earth3P-HR. The HighResMIP protocol (Haarsma et al., 2016) of having a control and a historical simulation starting from the same initial conditions was designed to minimize the consequences of such trends. Under the assumption that the model trend is similar for both simulations, it can be eliminated by subtracting the control from the historical simulation. Indeed, the global annual mean SAT and the AMOC of hist-1950 minus control-1950 display a very similar behavior in EC-Earth3P and EC-Earth3P-HR (Figs. 18b and 17c), with an upward trend for SAT and a downward trend for the AMOC. For SAT, the upward trend compares well with ERA5.

Figure 16Seasonal distribution of SSWs per decade in a [−10, 10] d window around the SSW date for ERA-Interim (black), EC-Earth3P (red) and EC-Earth3P-HR (blue) from control-1950. Time series are smoothed with an 11 d running mean. The total decadal frequency of SSWs is indicated in brackets.

Download

Weather regimes

Another way to test the representation of the midlatitude atmospheric flow, with a focus on the low-frequency variability (5–30 d), is to assess how well the models reproduce the winter (DJF) Euro-Atlantic weather regimes (Corti et al., 1999; Dawson et al., 2012).

Figure 17Time series of the annual AMOC index for the control-1950 (a) and hist-1950 (b) runs. Solid lines display the ensemble mean for the EC-Earth3P (red) and EC-Earth3P-HR (black). Shaded areas represent the dispersion due to the ensemble members. Black stars in panel (b) display values of RAPID data. (c) Mean ensemble difference between hist-1950 and control-1950 for Earth3P (red) and EC-Earth3P-HR (black).

Download

Figure 18Global mean averaged annual SAT (^∘C) in hist-1950 (a) for the three members of EC-Earth3P (red colors) and EC-Earth3P-HR (grey colors). (b) Mean ensemble difference between hist-1950 and control-1950 for EC-Earth3P (red) and EC-Earth3P-HR (black). ERA5 is indicated by the green curves. Panel (b) is scaled so that the starting point fits with the EC-Earth curves.

Download

Figure 19Observed cluster patterns for ERA (a), simulated cluster patterns in hist-1950 for EC-Earth3P (b) and EC-Earth3P-HR (c). The frequency of occurrence of each regime is shown above each subplot.

The analysis has been applied here to the EC-Earth3P and EC-Earth3P-HR hist-1950 simulations. Following recent works (Dawson and Palmer, 2015; Strommen et al., 2019), we computed the regimes via k-means clustering of daily geopotential height anomalies at 500 hPa over 30–85^∘ N, 80^∘ W–40^∘ E. As a reference, we considered the ECMWF reanalysis data from ERA40 (1957–1978) and ERA-Interim (1979–2014). The clustering is performed in the space spanned by the first four principal components obtained from the reference data set. More details on the technique used and on the metrics discussed here can be found in Fabiano et al. (2020) and references therein. Each row in Fig. 19 shows the resulting mean patterns of the four standard regimes – NAO+, Scandinavian blocking, Atlantic Ridge and NAO− – for ERA (top), EC-Earth3P (middle) and EC-Earth3P-HR (bottom). The regimes are quite well represented in both configurations. However, the matching is better in the standard-resolution version both in terms of rms and pattern correlation averaged over all regimes (see Table 3). Only the Scandinavian (Sc) blocking pattern is improved in the HR, whereas the other patterns are degraded. The most significant degradation is seen for the NAO− pattern, which is shifted westward in the HR simulation. The result for EC-Earth3P(-HR) goes in the opposite direction of what has been observed in Fabiano et al. (2020), where most models showed a tendency for improving the regime patterns with increased resolution. Concerning the regime frequencies, both model versions show a tendency to produce less NAO+ cases than the observations and more Atlantic Ridge cases (Fig. 19).

Table 3Some metrics to assess the overall performance in hist-1950 of the EC-Earth3P and EC-Earth3P-HR simulations in terms of weather regimes. The table shows the average root mean square error (RMSE) from the observed patterns and the relative average pattern correlation among all regimes, the average persistence of the two NAO states in days, the optimal ratio and the sharpness. The errors refer to the spread between members (standard deviation).

Download Print Version | Download XLSX

Another quantity of interest is the persistence of the regimes, since models usually are not able to reach the observed persistence of the NAO+/− states (Fabiano et al., 2020). As stated in Table 3, this is also observed for the EC-Earth3P hist-1950 simulations, and the effect of the HR is to increase the persistence of NAO+ but decrease that of NAO−.

Even if the HR is degrading the regime patterns, it produces a small but positive effect on the geometrical structure of the regimes. This is shown by the last two quantities in Table 3: the optimal ratio and the sharpness. The optimal ratio is the ratio between the mean inter-cluster squared distance and the mean intra-cluster variance: the larger the optimal ratio, the more clustered the data. The sharpness is an indicator of the statistical significance of the regime structure in the data set in comparison with a randomly sampled multi-normal distribution (Straus et al., 2007). The closer the value is to 100, the more significant the multimodality of the distribution. The sharpness tends to saturate at 100 for very long simulations, so the values reported in Table 3 are obtained from a bootstrap on 30 years chosen randomly. Both the optimal ratio and the sharpness are too low in the EC-Earth3P simulations, as is usually seen for all models. A significant increase with EC-Earth3P-HR is seen for the optimal ratio, and a smaller (non-significant) one is seen for the sharpness.

The increased resolution simulations have a stronger regime structure and are closer to the observations in this sense. However, the regime patterns are degraded in the HR version and this affects the resulting atmospheric flow. A similar result was obtained by Strommen et al. (2019) for a different version of EC-Earth and two other climate models.