# Description and Evaluation of the CNRM-Cerfacs Climate Prediction System (C3PS)

E. Sanchez-Gomez, Roland Séférian, Lauriane Batté, Sarah Berthet, Christophe Cassou, B. Dewitte, M. P. Moine M, R. Msadek, C. Prodhomme, Yeray Santana-Falcón, et al.

► To cite this version:

E. Sanchez-Gomez, Roland Séférian, Lauriane Batté, Sarah Berthet, Christophe Cassou, et al. De-scription and Evaluation of the CNRM-Cerfacs Climate Prediction System (C3PS). Journal of Ad-vances in Modeling Earth Systems, 2024, 16, 10.1029/2023MS004193. insu-04803119

HAL Id: insu-04803119

https://insu.hal.science/insu-04803119v1

Submitted on 25 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

10.1029/2023MS004193

**Key Points:**

* This study introduces and assesses C3PS, the new CNRM-Cerfacs Earth System-based prediction platform

* The platform can provide climate predictions from seasonal to multiannual timescales for relevant physical and biogeochemical fields

* The most outstanding result is the ability of C3PS to predict the net primary production and carbon fluxes at multiannual timescales

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

E. Sanchez-Gomez,
emilia.sanchez-gomez@meteo.fr

**Citation:**

Sanchez-Gomez, E., Séférían, R., Batté, L., Berthet, S., Cassou, C., Dewitte, B., et al. (2024). Description and evaluation of the CNRM-cerfacs climate prediction system (C3PS). *Journal of Advances in Modeling Earth Systems, 16*, e2023MS004193. https://doi.org/10.1029/2023MS004193

Received 22 DEC 2023
Accepted 3 OCT 2024

**Author Contributions:**

**Conceptualization:** E. Sanchez-Gomez, R. Séférían, L. Batté, S. Berthet, C. Cassou, R. Msadek, C. Prodhomme, A. Voldoire

**Formal analysis:** E. Sanchez-Gomez, R. Séférían, B. Dewitte, C. Prodhomme, Y. Santana-Falcón, L. Terray

**Funding acquisition:** E. Sanchez-Gomez, R. Séférían, L. Batté

**Investigation:** Y. Santana-Falcón

**Methodology:** E. Sanchez-Gomez, R. Séférían, L. Batté, S. Berthet, C. Cassou, M. P. Moine M, R. Msadek, C. Prodhomme, A. Voldoire

© 2024 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_4.png" /></figure>

# Description and Evaluation of the CNRM-Cerfacs Climate Prediction System (C3PS)

E. Sanchez-Gomez¹,² ID, R. Séférían³ ID, L. Batté² ID, S. Berthet³ ID, C. Cassou¹, B. Dewitte¹,⁴,⁵ ID, M. P. Moine M¹, R. Msadek¹, C. Prodhomme⁶, Y. Santana-Falcón³ ID, L. Terray¹ ID, and A. Voldoire² ID

¹CECI, Université de Toulouse, CNRS, Toulouse, France, ²Météo-France, Toulouse, France, ³CNRM, Université de Toulouse, Toulouse, France, ⁴Center of Advanced Studies in Arid Zones (CEAZA), Coquimbo, Chile, ⁵Departamento de Biología Marina, Facultad de Ciencias del Mar, Center for Ecology and Sustainable Management of Oceanic Islands (ESMOI), Universidad Católica del Norte, Antofagasta, Chile, ⁶Earth Scan, London, UK

**Abstract** The CNRM-Cerfacs Climate Prediction System (C3PS) is a new research modeling tool for performing climate reanalyzes and seasonal-to-multiannual predictions for a wide array of Earth system variables. C3PS is based on the CNRM-ESM2-1 model including interactive aerosols and stratospheric chemistry schemes as well as terrestrial and marine biogeochemistry enabling a comprehensive representation of the global carbon cycle. C3PS operates through a seamless coupled initialization for the atmosphere, land, ocean, sea ice and biogeochemistry components that allows a continuum of predictions across seasonal to multiannual time-scales. C3PS has also contributed to the Decadal Climate Prediction Project (DCPP-A) as part of the sixth Coupled Model Intercomparison Project (CMIP6). Here we describe the main characteristics of this novel Earth system-based prediction platform, including the methodological steps for obtaining initial states to produce forecasts. We evaluate the entire C3PS initialization procedure with the most up-to-date observations and reanalyzes over 1960–2021, and we discuss the overall performance of the system in the light of the lessons learned from previous and actual prediction platforms. Regarding the forecast skill, C3PS exhibits comparable seasonal predictive skill to other systems. At the multiannual scale, C3PS shows significant predictive skill in surface temperature during the first 2 years after initialization in several regions of the world. C3PS also exhibits potential predictive skill in Net primary production (NPP) and carbon fluxes several years in advance. This expands the possibility of applications of forecasting systems, such as the possibility of performing multiannual predictions of marine ecosystems and carbon cycle.

**Plain Language Summary** The study introduces and assesses the new climate prediction platform C3PS developed by the CNRM-Cerfacs modeling group in the framework of the H2020 TRIATLAS project. This prediction system is based on the latest version of the CNRM Earth system model, CNRM-ESM2.1, and was designed to produce predictions from seasonal to multiannual scales. C3PS is the result of the joint long-term effort of experts in seasonal and decadal forecasting and modellers of ocean physics and biogeochemistry within the CNRM-Cerfacs research group. An innovative aspect of our study is that it focuses on validating the initialization procedure, which is not often done in other studies presenting forecasting systems. We believe that the study of the reconstructions created to initialize the climate prediction systems is relevant, and even more so in the context of the new applications offered in the prediction of marine biogeochemistry and carbon fluxes. Regarding forecast skill, C3PS exhibits comparable seasonal predictive skill to other systems. On a multi-year scale, C3PS shows potential skill not only in physical climate variables, but also in NPP and carbon fluxes up to 3 years in advance, which extends the possibilities of application to marine ecosystems and multi-year carbon cycle forecasts.

## 1. Introduction

The field of near-term climate prediction has grown rapidly since the pioneering studies of Smith et al. (2007), Keenlyside et al. (2008), Pohlmann et al. (2009) and the very first attempt of decadal prediction coordinated experiments as conducted under the umbrella of the Fifth Phase of the Coupled Model Intercomparison Project (CMIP5). The analysis of CMIP5 decadal prediction experiments revealed a wide range of skill for different variables and across various prediction systems (Bellucci et al., 2015; Doblas-Reyes et al., 2013; Garcia-
**Software:** E. Sanchez-Gomez, R. Séférían, L. Batté, M. P. Moine M, A. Voldoire

**Validation:** E. Sanchez-Gomez, R. Séférían, L. Batté, S. Berthet, C. Cassou, B. Dewitte, R. Msadek, Y. Santana-Falcón

**Writing – original draft:** E. Sanchez-Gomez, R. Séférían, B. Dewitte, Y. Santana-Falcón, L. Terray

**Writing – review & editing:** E. Sanchez-Gomez, R. Séférían, L. Batté, C. Cassou, B. Dewitte, M. P. Moine M, R. Msadek, C. Prodhomme, Y. Santana-Falcón, L. Terray, A. Voldoire

Serrano et al., 2015 amongst others). Following this, CMIP6 undertook a new decadal prediction coordinated exercise with improvements with respect to CMIP5 (Boer et al., 2016). These improvements include not only model improvements, but also the increase of the number of starting dates and ensemble members in the decadal forecast archive, in order to ensure a robust assessment of decadal predictive skill. Results from CMIP6 show a substantial improvement of the Sea Surface Temperature (SST) prediction in the North Atlantic, in particular over the subpolar gyre (Borchert et al., 2021; Delgado-Torres et al., 2022). Over land, a significant increase in prediction skill of surface air temperature (SAT) is also reported (Monerie et al., 2018; Smith et al., 2019; Wu et al., 2019). Moreover, by the use of large ensembles, skillful predictions have been achieved for atmospheric patterns such as blocking (Athanasiadis et al., 2020; Schuster et al., 2019) and the North Atlantic Oscillation (Smith et al., 2020).

While the potential for useful applications has been demonstrated, the CMIP5/CMIP6 experiments have also highlighted a number of outstanding research questions and challenges in the climate prediction field (Bojovic et al., 2019; Cassou et al., 2018; Dunstone et al., 2022; Kirtman et al., 2013; Meehl et al., 2014; O’Kane et al., 2023). Previous decadal prediction exercises highlight the need for a better understanding of three key aspects for better exploiting the climate predictive potential and improving estimates of climate predictability at different timescales (Cassou et al., 2018; Keenlyside and Ba, 2010; Verfaillie et al., 2020): (a) The physical mechanisms of climate predictability, (b) initialization, and response to external forcing; (c) and an improvement of the forecast quality evaluation process.

One of the outstanding challenges is to identify the extent to which model prediction skill across a continuum of time-scales may benefit from initialization. Indeed, by establishing a framework for testing the added value of model initialization, as well as prescribing external forcings, decadal prediction systems have bridged the gap between well-established seasonal prediction and near-term projections (Meehl et al., 2009). In this sense, decadal predictions can provide seamless climate information from 1 month to several years ahead, offering the opportunity of exploring predictability at different timescales (Choi et Sun, 2022). This is relevant as it provides climate information addressing a growing demand from policy makers and stakeholders in the context of climate risk management.

Moreover, the required reduction of human-induced CO2 emissions and the need for adaptation of several sectors have, over the recent years, widened the range of application of climate predictions, with the inclusion of new Earth System components. Earth System Models (ESMs) have been recently implemented in climate prediction systems, allowing to explore the predictability of marine biogeochemistry and marine ecosystems (Park et al., 2019; Séférían et al., 2014; Yeager et al., 2022), terrestrial carbon fluxes (Séférían et al., 2018) and air-sea carbon fluxes and carbon budgets (Ilyna et al., 2021; Lovenduski et al., 2019). Following this path and in order to provide seamless seasonal to multiannual predictions for relevant physical and Earth system variables, the CNRM-Cerfacs modeling group has developed a new prediction platform, called C3PS, which is based on the CNRM-ESM2.1 model (Séférían et al., 2019). The birth of C3PS was possible by bringing together the expertise of the CNRM-Cerfacs modeling group in terms of seasonal and multiannual climate predictions and the latest developments in Earth system modeling made in the context of CMIP6 (Séférían et al., 2020).

In the present study, we introduce the C3PS system, by highlighting its main characteristics, and the seamless/coupled initialization method used for the atmosphere, ocean and marine biogeochemistry components. This coupled initialization has been achieved to enable the investigation of predictability across a continuum of time-scales, from seasons to years. Conversely to other studies presenting climate prediction systems, we perform an exhaustive evaluation of the initialization procedure, assessing its strengths and weaknesses. Finally, we also evaluate the performance of the C3PS system based on a variety of diagnostics and metrics.

Section 2 describes the main characteristics of C3PS, the initialization procedure and the experimental protocol used to perform the seasonal to multiannual predictions. Section 3 presents the reference data sets and metrics used. Section 4 provides a basic evaluation of the assimilation experiments used in the C3PS initialization. Section 5 assesses the skill of essential physical and biogeochemical fields at different time scales, and the concluding remarks are presented in Section 6.
## 2. Earth System-Based Prediction Platform

### 2.1. Model Description

The backbone of the C3PS platform is CNRM-ESM2-1 which is the Earth System model of second generation developed by CNRM-Cerfacs modeling group for CMIP6 (Séférrian et al., 2019).

The atmosphere component of CNRM-ESM2-1 is based on the global spectral model ARPEGE-Climat version 6.3 (Roehrig et al., 2020). ARPEGE-Climat resolves atmospheric dynamics using a T127 linear truncation. The physics is resolved on the corresponding reduced grid which offers a spatial resolution of about 150 km in both longitude and latitude. CNRM-ESM2-1 employs a ‘high-top’ configuration with 91 vertical levels that extend from the surface to 0.01 hPa in the mesosphere; 15 hybrid σ-pressure levels are available below 1,500 m.

The atmospheric chemistry scheme of CNRM-ESM2-1 is Reactive Processes Ruling the Ozone Budget in the Stratosphere version 2 (REPROBUS-C_v2). This scheme resolves the spatial distribution of 63 chemistry species but does not represent the low troposphere ozone non-methane hydrocarbon chemistry. CNRM-ESM2-1 also activates an interactive tropospheric aerosol scheme included in the atmospheric component ARPEGE-Climat. This aerosol scheme, named Tropospheric Aerosols for ClimaTe In CNRM (TACTIC_v2), represents the main anthropogenic and natural aerosol species of the troposphere.

The surface state variables and fluxes at the surface-atmosphere interface are simulated by the SURFEX modeling platform version 8.0 over the same grid and with the same time-step as the atmosphere model. Over the land surface, CNRM-ESM2-1 uses the ISBA-CTRIP land surface modeling system to solve energy, carbon and water budgets at the land surface (Decharme et al., 2019; Delire et al., 2020). Its physical core explicitly solves the one-dimensional Fourier and Darcy laws throughout the soil, accounting for the hydraulic and thermal properties of soil organic carbon. It uses a 12-layer snow model of intermediate complexity that allows to separate water and energy budgets for the soil and the snowpack. CTRIP is a dynamic river flooding scheme in which floodplains interact with the soil and the atmosphere through free-water evaporation, infiltration and precipitation interception. The ISBA-CTRIP land surface scheme also embeds a two-dimensional diffusive groundwater scheme to represent unconfined aquifers and upward capillarity fluxes into the superficial soil. More details on these physical aspects can be found in (Decharme et al., 2019). ISBA-CTRIP captures the land carbon cycle and vegetation-climate interactions with the representation of plant physiology, carbon allocation and turnover, and carbon cycling through litter and soil. It includes a module for wildfires, land use and land cover changes, and carbon leaching through the soil and transport of dissolved organic carbon to the ocean. A detailed description of the terrestrial carbon cycle can be found in Delire et al. (2019).

The ocean component of CNRM-ESM2-1 is the Nucleus for European Models of the Ocean (NEMO) version 3.6 (Madec et al., 2017) which is coupled to both the Global Experimental Leads and ice for ATmosphere and Ocean (GELATO) sea-ice model (Salas Mélia, 2002) version 6. NEMOv3.6 operates on the eORCA1L75 grid (Mathiot et al., 2017) which offers a nominal resolution of 1° to which a latitudinal grid refinement of 1/3° is added in the tropics; this grid describes 75 ocean vertical layers using a vertical z*-coordinate with partial step bathymetry formulation (Barnier et al., 2006).

The ocean biogeochemical component of CNRM-ESM2-1 uses the Pelagic Interaction Scheme for Carbon and Ecosystem Studies model version 2 coupled with trace gases module (PISCESv2-gas), which derives from PISCESv2 as described in Aumont et al. (2015). PISCESv2-gas simulates the distribution of five nutrients (from macronutrients: nitrate, ammonium, phosphate, and silicate to micronutrient: Iron), which regulate the growth of two explicit phytoplankton classes (nanophytoplankton and diatoms). PISCESv2-gas also simulates the ocean carbon cycle with the ocean carbonate chemistry, that is the dissolved inorganic carbon (DIC) and the alkalinity (Alk) and two organic carbon pools. The dissolved oxygen is prognostically simulated using two different oxygen-to-carbon ratios, one when ammonium is converted to or mineralized from organic matter, the other when oxygen is consumed during nitrification. Their values have been set respectively to 131/122 and 32/122. At the ocean surface, PISCESv2-gas exchanges carbon, oxygen, dimethylsulfide (DMS) and nitrous oxide (N₂O) tracers with the atmosphere using the revised air-sea exchange bulk formulation as in Wanninkhof (2014). PISCESv2-gas uses several boundary conditions which represent the supply of nutrients from five different sources: atmospheric deposition, rivers, sediment mobilization, sea-ice and hydrothermal vents.
## 2.2. Forcings

This section details the CMIP6 external forcing implementation into the C3PS platform. CNRM-ESM2.1 is run in concentration-driven mode to align with the CMIP6 Decadal Prediction Project (DCPP) protocol (Boer et al., 2016). For all the experiments whose simulated period lies within the historical period as labeled by CMIP6, that is, from 1850 to 2014, we implement the exact set-up that was used for the contribution to the CMIP6/DECK historical experiment (Eyring et al., 2016). Greenhouse gas concentrations, with the exception of stratospheric ozone, are implemented as recommended by Meinshausen et al. (2017). Specifically, in concentration-driven mode, atmospheric CO2 is restored to a global mean concentration and is not accounted for by the land and ocean carbon cycle modules. For details on the implementation of the forcings for CMIP6, please refer to Séférían et al. (2019) and Michou et al. (2020).

For simulated years after 2014 and in accordance with the DCPP protocol, the Shared Socioeconomic Pathway (SSP) 2–4.5 scenario forcing is prescribed (O’Neill et al., 2016). This is the “middle-of-the-road” scenario of the SSP2 socioeconomic pathways, with an intermediate 4.5 W/m² radiative forcing level by 2,100 (Gidden et al., 2019).

The major difference between the implementation of the external forcing in the C3PS platform and the usual CMIP6 simulation set-up for CNRM-ESM2-1 is the volcanic forcing. The CMIP6 experimental protocol now requires the use of a stratospheric volcanic background forcing (monthly climatology computed from years 1850 to 2000 volcanic forcing) during pre-industrial and future eras. However, over the 1850–2014 period, the volcanic forcing can be lower than the background forcing as used for the future period (beyond 2015). In consequence, we applied a linear ramp-up from the 2014 level to the background level over the 2015–2025 period, as suggested in Gillett et al. (2016).

## 2.3. Workflow and Data Production

Currently, the C3PS platform provides for both, seasonal and multiannual timescales the variables requested in the DCPP/CMIP6 tables (Boer et al., 2016), which are those variables relevant for forecast evaluation against observational data sets. Besides, we have included additional biogeochemical and ocean physics variables that are necessary to force marine ecosystem models. Most of these variables are already requested by the FishMIP initiative (Tittensor et al., 2019). Higher frequency variables, such as daily ocean potential temperature and oxygen are also. Concerning the atmosphere, C3PS also provides daily low-level winds (~100 m) and solar radiation variables as requested for renewable energy applications. The last mentioned variables are not available on ESGF but are available upon request.

The C3PS platform follows the DCPP/CMIP6 experimental protocol with regard to the multiannual predictions, although additional members have been performed to increase the ensemble size from 10 to 15 members.

All the C3PS related simulations were performed on the Belenos supercomputer, hosted at Météo-France site in Toulouse from June 2021 to February 2022. The work-flow is handled by the ECLIS (Environment for CLImate Simulations) package tool that was developed by the CNRM (https://www.umr-cnrm.fr/cm/spip.php?article14).

ECLIS is an ensemble of scripts and tools that allow for setting up and running all the experimental protocols performed by the CNRM-Cerfacs modeling group within CMIP and beyond. In particular, C3PS required additional ECLIS developments, such as dedicated scripts for the perturbation of initial atmospheric conditions and the management and launching of the members for all the starting dates (see Section 2.4).

The C3PS diagnostics production is managed by the XIOS output server (Meurdesoif, 2018; Meurdesoif, 2018). XIOS has been implemented in all the models developed by the CNRM-Cerfacs group, in particular to facilitate the huge CMIP6 data production. XIOS allows for declaring a priori the requested variables to be saved in the output files for a given experiment. Moreover, XIOS performs online operations on fields, such as spatial and vertical interpolations, vertical, spatial and time averages, vertical level extraction, thus saving a lot of post-processing time. XIOS has also been adapted to produce netCDF “CMOR” (Climate Model Output Rewriter) format files compliant with the CMIP6 Data Request specificities. More information about XIOS functioning can be found in Voldoire et al. (2019).
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematics of the initialization procedure of the C3PS platform. The “spring” connecting the two NEMO-PISCES boxes indicates that the coupled simulation is relaxed to potential temperature and salinity in the forced simulation.</figcaption></figure>

## 2.4. Seamless Prediction Procedure and Simulations

Most of the efforts involved in the development of the C3PS platform were oriented to achieve a satisfactory initialization procedure. In this regard, several challenges needed to be tackled. The first challenge was to participate in the DCPP/CMIP6, for which the required hindcast period starts in 1960, when biogeochemical observations required to initialize the biogeochemistry model are practically non-existent. A second challenge is how to robustly initialize a seamless climate prediction platform in which a continuum of timescales need to be considered. For seasonal prediction, atmospheric initialization is relevant for climate prediction (Materia et al., 2014). For longer timescales, atmospheric initialization is less relevant as the predictability mostly lies on the ocean and sea ice persistence and memory. A third challenge is to minimize the climate drifts that occur when the model is initialized from a state away from the climate model attractor. Besides, physical coherence amongst the initial states of all the model components of CNRM-ESM2-1 is necessary in order to avoid incompatibilities that could lead to abrupt initial shocks right after the initialization (Bilbao et al., 2021; Pohlmann et al., 2017; Sanchez-Gomez et al., 2016). Although model drifts in climate prediction systems are partly corrected before skill assessment, it is preferable to minimize them as much as possible to better distinguish the predictable signals (Meehl et al., 2022). For the initialization methodology, we built open the work of Sanchez-Gomez et al. (2016) which developed initial conditions for decadal predictions in CMIP5. This study provides a detailed analysis of the initial drifts and shocks experienced by a climate model when initialized for climate predictions. They also introduced an initialization methodology that involves nudging the ocean component of the coupled model toward reference data (see below). This approach was shown to reduce initial shocks and drifts and has been successfully applied in subsequent studies, including the decadal forecasting system based on the EC-Earth model in CMIP6 (Bilbao et al., 2021).

In order to overcome the three main challenges mentioned above, in the development of C3PS we have implemented an experimental protocol which is carried out in three main steps (Figure 1).

* Step 1: The initial idea was to use ORAS5 reanalysis (C3S 2021) as reference data for the ocean; however, a problem was identified in the first version of ORAS5 related to the non-stationarity of the Atlantic Meridional Overturning Circulation (AMOC), which degraded the predictive skill of the seasonal forecasts (Tietsche et al., 2020). Although a new version of ORAS5 is now available, it was not feasible to rerun all the simulations within the project's timeline. For this reasons, pseudo-observations are obtained through an ocean forced simulation in which the NEMO-PISCESv2gas model is forced by atmospheric fields from the JRA55do reanalysis (Tsujino et al., 2020) over the period 1960–2021 (Figure 1). This simulation (referred to as FORCED hereinafter) has been performed under the framework of the Global Carbon Project (GCP) (Friedlingstein et al., 2022; Hauck et al., 2020). A sea surface salinity relaxation is applied only beneath the sea ice, targeting the JRA55do climatology for the period 1955–2012. This relaxation uses a restoring coefficient of −166.67 mmd⁻¹, which is a standard value in the NEMO namelist for the forced configuration. The FORCED experiment was launched after a spin-up of 300 years in which the NEMO-PISCESv2gas model was forced by repeated cycles of 5 years corresponding to the 1958–1962 period. The analysis of this spin-up reveals that surface physical fields such as SST and salinity (SSS), and integrated fields such as ocean heat content (OHC) and AMOC are almost stabilized after the spin-up.
* Step 2: the 3D potential temperature and salinity fields issued from the FORCED simulation were used to constrain the ocean component of CNRM-ESM2.1 through a Newtonian damping procedure (Figure 1). This nudging simulation is performed over the period 1960–2021, and serves to generate the so-called dcppA-assim experiment according to the DCPP/CMIP6 experiment-id (Boer et al., 2016). The dcppA-assim (referred to ASSIM hereafter) can be considered as an in-house zero-order reanalysis product from which the initial conditions for all the components of CNRM-ESM2.1 are issued. As mentioned above, the methodology of the nudging was previously implemented and used in Sanchez-Gomez et al. (2016). It was shown to be beneficial to produce initial states physically consistent amongst all the components of CNRM-CM5, (b) to get initial states for the components with non-available observations and (c) to minimize the initial shock and drift in the prediction experiments. Here we use the same nudging strategy which consists in (a) a sea surface restoring of temperature and salinity of the NEMO component toward SST and SSS from the FORCED simulation; (b) a 3D Newtonian damping in temperature and salinity below the mixed layer to constraint the ocean subsurface toward FORCED. The sea surface restoring is applied globally in terms of heat and freshwater fluxes. The values of the restoring coefficients are −40 Wm⁻² K⁻¹ and −864 mmd⁻¹ for the heat and freshwater fluxes respectively. Note that the value of the coefficient for freshwater flux significantly differs for those used in previous studies (Bilbao et al., 2021; Servonnat et al., 2015; Sanchez-Gomez et al., 2016). The rationale of this is to have the same restoring time scale for SST and SSS, that is 60 days for a mixed layer of 50 m (Barnier et al., 1995). The 3D Newtonian damping is applied as follows: On the vertical, there is no damping above the mixed layer to allow for physical coherence between the mixed layer and the surface processes. Below the mixed layer down to 800 m depth, the damping term is set to 10 days and for the deep ocean below, a weak damping is used (~1 year). Horizontally, subsurface nudging is only applied outside the 15°S–15°N latitudinal band and from 300 km off the coast to avoid spurious vertical currents at the equator and coastal effects respectively (Sanchez-Gomez et al., 2016). A buffer zone of 5°, which decreases linearly, is applied between the nudged areas and the rest of the ocean. Similar nudging methodologies are also adopted in Bilbao et al. (2021) in order to obtain initial states for seasonal and decadal predictions. The ASSIM simulation has been duplicated with a set of perturbed parameters in order to obtain an ensemble of three members. For this, the ocean and atmosphere diffusivity have been slightly perturbed separately to produce additional ASSIM members (see Figure 1).

* Step 3: the ASSIM ensemble will be used as initial conditions for all the CNRM-ESM2.1 components for both seasonal and multiannual predictions. Only the atmospheric restarts provided by ASSIM are modified in order to adapt C3PS to seasonal forecasting. For this purpose, the dynamical fields contained in the restarts of ARPEGE in the ASSIM ensemble are replaced by the dynamical fields provided by the ERA5 reanalysis (Hersbach et al., 2020). Finally, the prediction procedure is performed as follows: For the seasonal timescale, two initializations per year are considered, that is, 1st May and 1st November. For each start date, an ensemble of 30 members is generated. The atmosphere is perturbed by using a small increment of the atmospheric dynamical fields provided by ERA5. This increment, introduced only at the initialization time, is drawn randomly from a set of increments computed during a previous historical atmospheric nudging simulation where the ARPEGE model is weakly constrained toward the ERA5 reanalysis (Batté & Déqué, 2012). Ten increments were used for each ASSIM member, thus building a 30-member ensemble. Seasonal predictions starting 1st May are run for 6 months. For the multiannual timescale, the perturbation procedure is identical to that of the seasonal scale, except that only the forecast starting on 1st November is continued up to 5 years, and with only 15 members. Hereinafter the set of seasonal to multiannual predictions will be referred to as PRED.

The most significant differences from the previous work by Sanchez-Gomez et al. (2016) in terms of the initialization methodology are: (a) The use of a forced ocean simulation with biogeochemistry instead of an ocean reanalysis, and (b) the initialization of the atmosphere, by incorporating ERA5 dynamical fields into the *e* atmospheric component restart.

# 3. Data Sets and Methods to Assess C3PS Performances

## 3.1. Reference Data Sets for Verification

Several observational and pseudo-observational products have been used to evaluate the ASSIM reconstruction (Section 4) and to compute the forecast skill scores for PRED (Section 5).
The physical variables we have considered are: air temperature at 2 m (SAT), ocean temperature and salinity, OHC, Arctic sea ice concentration (SIC) and extent (SIE) and AMOC. To evaluate ASSIM sea surface temperature, we use a blended product consisting of an average of the Hadley Center Sea Ice and SST version 1 (HadISST1, Rayner et al., 2003) and ERST v5 (Huang et al., 2017) over ice-free sea water. Over land and over sea-ice we average BEST (Muller, Curry, et al., 2013; Rohde et al., 2013), CRU-TS4-00 (Harris et al., 2014), and GHCN-CAMS (Fan and van den Dool, 2008).

The latest EN4 (version 4.2.2) objective analysis product is used as a reference for 3D ocean temperature and salinity (Good et al., 2013) and for OHC computation. This is a 1° × 1° gridded data set derived from ocean and temperature profiles with quality checks, which runs from 1900 to present. Here we have considered the EN4 analyses with the Gouretski and Reseghetti (2010) bias correction. SST and SIC reference data are also issued from HADISST1. We use the RAPID time series of the AMOC measured at 26°N as reference data (Moat et al., 2022), which are available from 2004 to 2020.

To analyze biogeochemistry, we focus on surface chlorophyll, integrated net primary production and global (land and ocean) carbon fluxes. Monthly means of chlorophyll-a concentration with a spatial resolution of 1° were issued from the ESA Ocean Color Climate Change Initiative (ESA-OC-CCIV3.1) project (Valente et al., 2022, https://climate.esa.int/en/projects/ocean-colour/). Net primary production (NPP) was obtained using a spectrally resolved model to simulate changes in photosynthesis as a function of irradiance (Kulk et al., 2020). This model incorporates vertical structure in chlorophyll-a concentration from OC-CCIV4.1. NPP data are 1° gridded and are available for the period 1998–2021. Carbon fluxes are evaluated by using the GCP reconstruction between 1959 and 2021 (Friedlingstein et al., 2022). This reconstruction currently represents the best estimates of the global carbon sink over the industrial era since 1959. For the ocean carbon sink (fCO₂) we use the Surface Ocean CO₂ Atlas version 2022 (SOCATv2022; Bakker et al., 2022) for the period 1990–2021.

To evaluate ASSIM reconstruction, besides the observational and analysis products described above, we consider the FORCED simulation, and the historical experiment performed with CNRM-ESM2.1 for CMIP6/DECK (Séférían et al., 2019) and referred here as FREE, which represents the free model (no data assimilation) run. The FREE ensemble consists of 10 members.

## 3.2. Metrics for Skill Assessment

The full-field initialization strategy used in C3PS requires to remove the forecast drift that inevitably occurs in any climate prediction system before performing the verification with observations and the skill estimate. We use the standard approach of transforming the raw model data into anomalies relative to the climatological forecast for each lead time,

$$X'_{j,l} = X_{j,l} - X_l \qquad (1)$$

where $X_{j,l}$ represents the ensemble-mean forecast from starting date j at lead time l and $X_l$ is the average over these forecasts over all starting dates for a given lead time. This is the so-called mean drift correction method, which assumes that forecast drift does not depend on the background climate state, that is, the drift is not considered to change between two different climate states from the point of view of global warming (Garcia Serrano & Doblas-Reyes, 2012; Meehl et al., 2014). Note that for the forecast period 1960–2021, the number of starting dates is 62 × 2 for the seasonal, and 62 for the multiannual timescales. In the case of multiannual forecasts, starting on 1st November each year, we focus our analysis on the following 5 years beginning in January (2 months after the initialization).

For both seasonal and multiannual timescales we use the standard verification framework as outlined in Goddard et al. (2013). We rely on the anomaly correlation coefficient (ACC), root mean square error (RMSE) and the Mean Square Skill Score (MSSS). The MSSS is especially used to assess the added value of the initialization and it is computed following equations 4–6 from Goddard et al. (2013). A MSSS score greater than 0 means that PRED is more accurate than FREE. For the seasonal forecast, persistence scores are used as a benchmark of C3PS scores and a t-test is used for assessing the statistical significance of the correlation.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Departure in blended surface temperature of FREE (a) and ASSIM (b) simulations from observations over 1960–2014. Blended surface temperature combines surface-air temperature over land and sea ice and sea surface temperature over ice-free sea water. Observations average several data sets: HadISST1 (Rayner, 2003) and ERSST v5 (Huang et al., 2017) over ice-free sea water; BEST (Muller, Curry, et al., 2013; Rohde et al., 2013), CRU-TS4-00 (Harris et al., 2014), and GHCN-CAMS (Fan & van den Dool, 2008) over land and sea ice. Units are in degrees Celsius.</figcaption></figure>

According to Goddard et al. (2013), for the skill maps and in order to remove small-scale unpredictable noise, all model and observational data are interpolated to a common 5-degree regular grid using the ESMF patch interpolation included in The Ncar Command Language (2019) –NCL.

We assess the added value of the initialization in C3PS by comparing the hindcasts PRED and the non-initialized historical ensemble (FREE) against FORCED or JRA55do for atmospheric variables. To properly evaluate skill differences between PRED and FREE, either through ACC or MSSS, a non-parametric bootstrap technique is used to assess the statistical significance of the skill scores (Goddard et al., 2013; Yeager et al., 2018). A block-bootstrap distribution of the scores is constructed at each location (grid point or time series) by resampling (with replacement) pairs of observations and hindcasts across the time dimension, and in addition, the PRED and FREE ensembles across the ensemble member dimension. Following these previous papers, we use a block size of 6 years as a trade-off between autocorrelation of the physical variables and the number of blocks (results are very similar to those based on 5 or 7-year blocks). The derived *p*-values are estimated as in Yeager et al. (2018).

Finally, the hindcast performance is evaluated by considering the so-called “potential predictability,” which as defined here consists of using as reference data set the FORCED and ASSIM experiments (Yeager et al., 2022). The skill calculated with respect to ASSIM represents the ability of the model to predict its own future data-constrained behavior. The notion of “potential predictability” is also interesting to assess forecast performance for biogeochemistry, since observations are available over a short period of time. Here, we will compare “potential predictability” versus the predictability with respect to FORCED for ocean and biogeochemistry, and with respect to JRA55do reanalysis for atmospheric variables (i.e., SAT).

## 4. Basic Evaluation of the C3PS Initialization Procedure

The assessment of the C3PS initialization strategy aims to determine how the nudging of the ocean physics has affected the performance of CNRM-ESM2.1 at simulating relevant physical and biogeochemistry fields.

Figure 2 shows that FREE exhibits common coupled model biases in the North Atlantic Ocean (the so-called “blue spot”), the southeastern Tropical Atlantic along the Benguela coast, and the equatorial Pacific cold tongue. Biases over the eastern ocean basins and Southern Ocean are reduced in ASSIM, as expected. Over land, ASSIM and FREE do not differ much in terms of biases, though over some regions like North America, Northern Africa temperature biases are slightly reduced. This fact indicates that ocean nudging does not have much impact over the continental areas.

Ocean temperature and salinity fields were used to nudge the ocean component of CNRM-ESM2.1 in order to generate ASSIM as explained above. Therefore, it is essential to evaluate both the performance of FORCED to simulate the mean state of the subsurface ocean as informed by observations, and then to evaluate the bias reduction of the ASSIM simulation with respect to the FREE simulation. Figure 3 displays thus the departure in ocean temperature and salinity at 100 m depth of FORCED with respect to observations, and of FREE and ASSIM with respect to FORCED. The FORCED simulation captures the main distribution of ocean temperature at the
<figure><img src="image_7.png" /><figcaption>**Figure 3.** (a, d) Difference of ocean temperature and salinity at 100 m depth between FORCED with respect to EN4 observations over 1960–2014. Differences between FREE (b, e) and ASSIM (c, f) with respect to FORCED over the same period. Spatial correlations and root mean square error (RMSE) of the time average over the whole period are shown on the top of each figure. Correlations and RMSE are computed against EN4 for (a, d) and against FORCED simulation for (b, c, e, and f). Observations are extracted from the quality-controlled EN4 data set (Good et al., 2013). Units are in degrees Celsius for temperature and psu for salinity. pCorr and pRMSE indicate that both correlations and RMSE are computed against the FREE simulation, rather than against observations, which is the case for Corr and RMSE.</figcaption></figure>

subsurface as depicted from observations. Nonetheless, FORCED overestimates temperature in the tropical Atlantic and across the North Pacific and the Southern Oceans, while it underestimates it southeast of New Zealand and in the tropical Pacific. In addition, the FORCED simulation strongly underestimates temperature in the North Atlantic where a well-documented “warming hole” has been related to a persistent slowdown of the AMOC (Drijfhout et al., 2012; Menary et al., 2018; Swingedow et al., 2021). By contrast, salinity is underestimated in this region, over the polar region and most regions of the Pacific Ocean, as reported in Voldoire et al. (2019).

As mentioned above, the nudging in subsurface waters is only applied in latitudes higher than ±15°. Accordingly, as seen in Figure 3, main differences of the ASSIM simulation with respect to the FORCED simulation occur in tropical regions, where the ASSIM tends to underestimate both temperature and salinity. In contrast, the underestimation of temperature expands to the Atlantic and South Pacific Oceans in the FREE simulation, while it overestimates temperature in the Southern Ocean and the California Current. The FREE simulation also underestimates salinity across most of the Atlantic and Pacific oceans, while it overestimates it across the North Pacific and the Indian oceans.

In conclusion, as expected, in both FORCED and ASSIM, the biases in surface and subsurface are strongly reduced compared to the FREE run, which confirms the validity of the methodology to generate oceanic initial conditions.

## 4.1. Drivers of Seasonal Climate Variability

In order to evaluate the realism of ASSIM in accounting for ENSO variability, we focus on how the nudging procedure impacts the ENSO diversity considering that this is a fundamental ENSO property that determines its seasonal evolution and teleconnections (Capotondi et al., 2020). The ENSO diversity or complexity (Timmerman et al., 2018) refers to the existence of warm and cold events with different SST patterns and amplitudes, with the extreme warm events being of Eastern Pacific type, while moderate warm and cold events being of Central Pacific type. Although the CNRM-ESM2.1 model (FREE) has some skill in simulating ENSO feedback strength (Lee et al., 2021), it has difficulty in simulating ENSO amplitude diversity, which manifests as a negative skewness of SST anomaly in the eastern equatorial Pacific.

Here as a compact measure of ENSO diversity (or nonlinearity), we use the value of the first coefficient of a quadratic fit in the phase plane of the first and second principal components (PCs) of SST anomalies in the tropical
<figure><img src="image_8.png" /></figure>

**Figure 4.** Phase space of the first and second principal components (PCs) (PC) of monthly Sea Surface Temperature anomalies in the tropical Pacific (120°E-80°W; 11°S-11°N) for observations (blue dots, from HadISST 1960-2020) and the ASSIM runs (orange and red dots). Nonlinearity is measured by fitted quadratic curves between PC values (blue: observations, red: ASSIM, black: FORCED, green: FREE). The PC axes have been rotated by 45° to infer the Eastern (E) and Central (C) Pacific indices. Three different types of observed El Niño events are highlighted with light blue circles (December): 1997 extreme Eastern Pacific El Niño, 2009 Central Pacific El Niño and 2015 mixed-type.

Pacific (Cai et al., 2018; Karamperidou et al., 2017), hereafter referred to as *α*. For HadISST1 data, the two branches of this quadratic fit tend to align along axis that correspond to the PC1 and PC2 axes rotated by 45° (Figure 4). The rotation of the PC time series defines the E and C indices, with $E = (PC1 - PC2)/\sqrt{2}$ and $C = (PC1 + PC2)/\sqrt{2}$, that account for the variability of Eastern Pacific events and Central Pacific events (Figure 5, top), respectively (Takahashi et al., 2011). While $\alpha = -0.33$ for observations, $\alpha = 0.10 \pm 0.06$ for FREE (the error corresponds to the standard deviation amongst the 10 members), which results from the negative ENSO asymmetry of the CNRM-ESM2.1 model (Lee et al., 2021). ASSIM has a more realistic ENSO non-linearity ($\alpha = -0.29$), almost identical to FORCED ($\alpha = -0.28$), indicating that the nudging procedure succeeds in restoring positive ENSO asymmetry to the observed value, along with improving ENSO diversity (see the blue curve paralleling the red curve in Figure 4). Still, ASSIM tends to have a larger ENSO variability than in the observations as evidenced by the larger amplitude of the E and C mode patterns compared to observations (Figure 5 middle), which is due to FORCED overestimating ENSO variability.

## 4.2. Drivers of Interannual to Decadal Climate Variability

To examine the drivers of interannual to decadal Pacific variability simulated by ASSIM, we focus on the Tripole Pacific Index (TPI) as defined by Henley et al. (2015). The TPI is a proxy of the Interdecadal Pacific Variability (IPV) and it is based on the difference between the SST anomalies averaged over the central equatorial Pacific minus the average of the SST anomalies in the

Northwest and Southwest Pacific (see Bilbao et al., 2021; Henley et al., 2015 for more details). Here, we do not consider SST anomalies as we are interested not only in the phase of the low frequency variability, but also in the model mean state. ASSIM and FORCED are coherent with HadISST1 SSTs evolution (Figure 6a), which is expected due to the sea surface restoring. The ensemble mean of the temporal correlation between each ASSIM member and FORCED with HadISST1 is 0.92 (see Table 1). Interannual variability of TPI is underestimated by the FREE ensemble as shown by Figure 6a and the variance ratio in Table 1. The smaller amplitude of Pacific decadal variability in the CNRM-Cerfacs models was also reported in Voldoire et al. (2019), which suggests a lack of the ENSO teleconnection at decadal timescales. In terms of RMSE, ASSIM presents an improvement with respect to FREE (Table 1). Interannual variability of OHC integrated over the first 300 m (OHC300) indicates that ASSIM is quite in phase with EN4, with a correlation value of 0.79 (Figure 6b and Table 1). ASSIM also improves the amplitude of the interannual variability with respect to FREE (see Table 1).

Regarding the OHC300 mean state, ASSIM and FORCED exhibit a cold bias which is weaker than in FREE (Table 1, RMSE). This cold bias of CNRM-ESM2.1 is also present in the coupled ocean-atmosphere climate model CNRM-CM6.1 (Voldoire et al., 2019). In general, like most coupled models, CNRM-Cerfacs models show a cold temperature bias in the Pacific Ocean from the surface to around 300 m depth. This cold bias is suggested to be caused by too strong surface winds curl exerting a pronounced wind stress curl into the surface ocean (see also Figure 3c). From Figure 6b, ASSIM mean state lies in between the reference data EN4 and FREE, indicating that initializing the ocean component of CNRM-ESM2.1 from ASSIM could potentially reduce the model drift in the predictions, which is actually the scope of our initialization procedure.

Another driver of interannual to decadal ocean variability is the Atlantic Multidecadal Variability (AMV). It was shown that CNRM-Cerfacs models simulate quite well the AMV spatial pattern in comparison to observations (Voldoire et al., 2019). Here we analyze the Subpolar Gyre in the North Atlantic (SPNA), which is closely correlated to the AMV. The SPNA SST time series (Figure 6c) exhibits a high temporal correlation with HadISST1 in both ASSIM and FORCED, with correlations of 0.90 and 0.93, respectively (see Table 1). Moreover, observations lie within the FREE multi-member spread, indicating that in terms of mean state, the free model performs quite well for this area. Note that the members of FREE show a pronounced variability, as also indicated in the variance ratio in Table 1. Indeed, the models CNRM-CM6.1 and CNRM-ESM2.1 are
<figure><img src="image_9.png" /><figcaption>**Figure 5.** C (left) and E (right) mode patterns for observations (top), FORCED and ASSIM (middle panels) and FREE (bottom). Dispersion (rms amongst the ensemble) is indicated for FREE in white contours.</figcaption></figure>

characterized by a large SST variance over the SPNA at decadal timescales, which is strongly correlated to AMOC variations, Arctic freshwater flux balance and northward salt transports from the tropical area (Voldoire et al., 2019). The marked decadal variability in FREE is also visible in the OHC integrated over the first 700 m (Figure 6d and Table 1). Once again, the correlation of ASSIM with regards to EN4 (0.91) indicates a good temporal coherency in the ocean subsurface.

As mentioned above, AMOC variations simulated by CNRM-ESM2.1 are highly correlated to decadal variability over the SPNA and Northern Seas (Voldoire et al., 2019). Time series of maximum AMOC at 26°N show a large low frequency variability in the members of FREE (Figure 7a), previously documented in Séférian et al. (2019) and Waldman et al. (2021). The mean AMOC value at 26°N of FREE is $16.4 \pm 2.3$ Sv for the period 1960–2014. The uncertainty in the latter value is estimated by considering one standard deviation amongst the members of the ensemble. The FREE ensemble AMOC is in good agreement with the RAPID mean value of 16.8 Sv for the observed period. Moreover, the depth of maximum observed AMOC is well simulated by FREE (Figure 7b). FORCED and ASSIM show a weaker AMOC (Figures 7a and 7b), with mean values of $12.2 \pm 0.7$ and $12.6 \pm 0.8$ respectively. The GCP experimental protocol used to perform FORCED is quite similar to that proposed in OMIP2/CMIP6 (Tsujino et al., 2020). The latter study documents that, in general, the forced ocean simulations show a lower AMOC intensity compared to RAPID. This underestimation of the AMOC is even more pronounced in the NEMO3.6/GELATO forced model configurations, suggesting that coupling with the atmosphere plays an important role in this high variability and intensity of the AMOC in the CNRM-Cerfacs models. The
<figure><img src="image_10.png" /><figcaption>**Figure 6.** (a) Tripole Pacific Index (TPI) annual time series from 1960 to 2021 for the Sea Surface Temperature (SST) and (b) ocean heat content (OHC) integrated over the first 300 m for the FREE ensemble (green), ASSIM ensemble (red), FORCED (black) and HadISST1/EN4 (blue). (c) Subpolar North Atlantic (SPNA) index annual time series from 1960 to 2021 for the SST and (d) OHC integrated over the first 700 m for the same experiments. For the OHC the observational reference is EN4. The TPI index is computed from raw data according to Henley et al. (2015). The SPNA index from raw data is obtained according to Bilbao et al. (2021) (SPNA: 50–65°N, 60–10°W). For FREE and ASSIM the ensemble means (thick line) and plus/minus one standard inter-members deviation is shown (red and green shading).</figcaption></figure>

nudging of temperature and salinity impacts the AMOC in ASSIM by altering density-driven deep water formation, leading to a correlation of 0.60 compared to FORCED.

The reason of the AMOC underestimation of ASSIM and FORCED can be partially explained by less dense subsurface waters of ASSIM and FORCED compared with FREE over the deep convection areas, that is, Labrador and GIN seas (Figure 7cd). These differences of density are mainly explained by warmer and less salty waters in FORCED and ASSIM, which are less realistic than those of FREE. The impact of the T/S nudging of CNRM-ESM2.1 toward FORCED seems to affect freshwater fluxes over the Labrador and GIN-Sea regions, since ASSIM is less salty than FORCED. The AMOC and related ocean deep convection characteristics in the ASSIM simulations are consistent with regional features of the Arctic SIC climatology (not shown). Indeed, FREE has a more extended sea ice area than ASSIM over the marginal seas in winter, which is consistent with colder and saltier waters over the Labrador and GIN Seas.

<table><caption>**Table 1**

Performance Metrics (Correlation, Variance Ratio and Root Mean Square Error) Computed With Respect to the Observational References and the Different Experiments: FORCED, ASSIM and FREE</caption><thead><tr><th></th><th>TPI (SST)<br>FORCED ASSIM FREE</th><th>TPI(OHC300)<br>FORCED ASSIM FREE</th><th>SPNA(SST)<br>FORCED ASSIM FREE</th><th>SPNA(OHC700)<br>FORCED ASSIM FREE</th></tr></thead><tbody><tr><th>Correlation</th><td>0.92 0.92 0.01</td><td>0.95 0.79 0.14</td><td>0.93 0.90 0.16</td><td>0.93 0.91 0.16</td></tr><tr><th>Variance ratio</th><td>1.17 1.06 0.83</td><td>0.96 0.83 0.71</td><td>0.91 0.84 1.30</td><td>0.95 0.90 1.57</td></tr><tr><th>RMSE</th><td>0.39 0.25 0.79</td><td>0.37 0.37 0.57</td><td>0.58 0.52 0.99</td><td>0.001 0.001 0.002</td></tr></tbody></table>

*Note.* The time series used to compute the metrics are displayed in Figure 6. The values shown for ASSIM and FREE are the ensemble mean of the values computed for each individual member.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** (a) Time series of the maximum Atlantic Meridional Overturning Circulation (AMOC) at 26°N for the FREE ensemble (green), ASSIM ensemble (red), FORCED (black) and RAPID data (blue). Units in Sv. (b) Vertical profile of AMOC at 26°N for the FREE ensemble (green), ASSIM ensemble (red), FORCED (black) and RAPID data (blue). Units in Sv. For (a and b) the FREE and ASSIM the ensemble means are shown (thick line) together with plus/minus one standard inter-members deviation (shading). (c) Temperature-Salinity diagram over the Labrador Sea area (70°W-45°W, 50°N-68°N) at 700 m depth for the FREE ensemble (green), ASSIM ensemble (red), FORCED (black) and EN4 data (blue). Panel (d) The same as panel (c) but for the GIN-Sea area (25°W-10°E, 65°N-80°N). Only ensemble means are shown for FREE and ASSIM. Potential density is computed from the NCL function "rho_mwjf."</figcaption></figure>

The annual cycle of SIE and SIV shows that ASSIM is comparable to the FREE ensemble (Figure 8ab), except from October to December where ASSIM performs better than FREE. In general, FORCED, ASSIM and FREE overestimate the maximum Arctic SIE which is connected to a too cold mean state with respect to HadISST1 (Figure 8a). The SIV simulated by ASSIM overlaps the FREE climatology, indicating a weak control of the

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Seasonal cycles of the Arctic Sea Ice Extent (SIE) (a) and Volume (b) computed in the period 1960 to 2021 for the FREE (green), ASSIM (red), FORCED (black) and HadISST1 (blue). For FREE and ASSIM the ensemble means (thick line) and plus/minus one standard inter-members deviation is shown (shading).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 9.** (a, d) Difference of ocean surface chlorophyll maximum (top panels) and minimum (bottom panels) between FORCED with respect to ESA-OC-CCIV3.1 observations over the period 1998 to 2017. Differences between FREE (b, e) and ASSIM (c, f) with respect to FORCED over the same period. Global average spatial correlations and root mean square error (RMSE) are shown on the top of each figure. Correlations and RMSE are computed against WOA2018 for (a, d) and against FORCED simulation for (b, c, e, and f). Surface chlorophyll maximum corresponds to the average over the months March, April, and May. Surface chlorophyll minimum corresponds to the average over the months August, September, and October. Observations correspond to monthly climatological data extracted from the quality-controlled 1° resolution ESA-OC-CCIV3.1 data set (Valente et al., 2022).</figcaption></figure>

nudging on the volume. The correlations between ASSIM and FORCED interannual time-series of SIE are 0.82 and 0.36 for March and September respectively (not shown). The strong correlation in March suggests that nudging significantly constrains SIE in the ASSIM ensemble during the winter freeze-up. In contrast, the weaker correlation in September indicates that this constraint is less effective during the melt season, which aligns with the greater influence of ocean conditions on sea ice during winter. Our results indicate that CNRM-ESM2.1 will be initialized from sea ice climatological conditions similar to those in the FREE experiment. This could be advantageous for modeling sea ice drift, which may otherwise negatively impact predictability in the SPNA region (Bilbao et al., 2021; Huang et al., 2015).

### 4.3. Biogeochemistry

The biases of both ocean surface chlorophyll maximum and minimum show that the FORCED simulation has difficulties in representing surface chlorophyll patterns (Figures 9a–9d). In general, the FORCED simulation underestimates the maximum chlorophyll values in the North Atlantic and North Pacific Oceans, while it overestimates both maximum and minimum chlorophyll observations in both the Pacific and Southern Oceans. The FORCED simulation also overestimates observations in the North Atlantic. The FREE biases with respect to the FORCED simulation are stronger over the western boundaries in the northern oceans and over the Southern Ocean for chlorophyll minimum (Figures 9b–9e). The difficulties of CNRM-ESM2.1 to represent surface chlorophyll over the Southern Ocean were documented in Séférian et al. (2019), and related to erroneous phytoplankton growth representation over the high-nutrients areas. Coastal chlorophyll biases were explained by deficiencies in remote-sensing products to represent coastal concentrations of surface chlorophyll (e.g., Gregg & Casey, 2004).

ASSIM biases with respect to the FORCED simulation are still similar to those shown for the FREE simulation (Figures 9c–9f). A lower RMSE and higher pattern correlation quantitatively indicate that ASSIM deviations from FORCED are smaller, pointing at a marginal impact of the sea surface restoring and 3D nudging. However, the nudging applied to temperature and salinity does not improve chlorophyll concentrations because it fails at improving distribution of nutrients in most regions. Indeed, an analysis on the biases of both surface nitrate (NO3) concentrations and mixed layer depth (MLD) between FREE and ASSIM with respect to WOA2018 climatology (Figure S1 in Supporting Information S1), suggests that the biases in NO3 are too strong to be compensated by the nudging on ocean physics. Moreover, in both the Southern Ocean and the North Atlantic, an underestimation of
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Scatterplot of the variance ratio experiment/observation versus temporal correlation (experiment/observations) for the integrated primary productivity averaged over the (a) whole tropical regions (30S-30°N), (b) the tropical Pacific, (c) the tropical Atlantic and (d) the Indian Ocean for the FREE ensemble (green), ASSIM ensemble (red), FORCED (black). Blue dashed line indicates a perfect match (=1) for the variance ratio. The observational reference is issued from Kulk et al. (2020) data set for the available period 1998–2018.</figcaption></figure>

the MLD together with an overestimation of NO3 may explain the consistent overestimation of surface chlorophyll minimum in those regions. Nutrient-rich waters that concentrate within a shallow MLD will strengthen the excessive development of phytoplankton. Phytoplankton growth in these regions will become limited by light availability, which explains why the overestimation of surface chlorophyll maximum is not as high as for surface chlorophyll minimum, especially over the North Atlantic.

The impact of the nudging on NPP is diagnosed over the tropical oceans in terms of interannual variability and temporal coherence with respect to observational estimates. Figure 10a shows that the nudging leads to an improvement of simulated interannual variance and temporal coherence with observations for ASSIM with respect to FREE in the tropical band. This improvement comes from the Tropical Pacific (Figure 10b), for which the correlation between ASSIM and observations is high (around 0.6). The fact that the sea surface restoring improves the phasing of the NPP interannual variability of CNRM-ESM2.1 with respect to NPP observational estimates over the Pacific Ocean was also documented by Séférrian et al. (2014). The SST restoring induces an improvement of SST gradients and short-term dynamical adjustment of winds, which combined with a good representation of nutrients over the area by CNRM-ESM2.1 can lead to a better simulated NPP. Contrary to the Pacific, oceanic nudging does not induce a clear impact in the NPP representation in the other tropical basins. In the Atlantic and the Indian oceans, ASSIM and FREE results are very similar.

At this point, the quality of FORCED and ASSIM for the biogeochemistry against observations may be questionable. Further evaluation of these reconstructions could be considered, but validation would be challenging due to the uncertainties and limitations of the observational biogeochemical data. While studying the variability of these marine biogeochemical variables is scientifically relevant, we believe it is beyond the scope of this paper, which focuses on introducing a novel prediction platform and documenting its performance. The performance of global biogeochemical models in ocean-forced configurations has been examined in several other publications stemming from the GCP/RECCAP2 initiative (Doney et al., 2024; Terhars et al., 2024). Although not without imperfections, the ocean and biogeochemical forced configurations used as surrogates for observations are considered state-of-the-art and are well-documented and assessed in the scientific literature.
<figure><img src="image_15.png" /><figcaption>**Figure 11.** (a) Annual time-series of the ocean carbon sink from 1990 to 2021 for FREE and ASSIM ensembles and Global Carbon Project (GCP) data product (Friedlingstein et al., 2022). The ocean carbon sink is represented in anomaly with respect to the long-term mean over the 1990–2021 period. The ensemble mean of available GCP ocean biogeochemical models and observational data products are given in gray and dark blue. For the sake of discussion, the ensemble of the eight available data-products is split in two sub-ensemble characterized by either a stronger (GCP data Strong Var, +) and lower (GCP data Low Var, ~) variability than the ensemble mean. (b) Scatter plot comparing model properties in terms of variability of the ocean carbon sink (y-axis) and the chronology of the ocean CO2 fugacity (fCO2) over the 1990–2021 period is provided for individual realization of FREE (green), ASSIM (red) and GCP models (gray). The ensemble average is given by the green, red and black crosses for FREE, ASSIM and GCP models. The dimensions of the vertical and horizontal bars on the crosses represent the magnitude of the spread of the GCP, ASSIM, and FREE ensembles for each axis.</figcaption></figure>

The impact of the nudging on the ocean carbon sink is assessed in terms of trends and variability in Figure 11. Figure 11a shows that FREE and ASSIM simulations capture the long-term increase of the global carbon sink as shown by the GCP reconstruction between 1990 and 2021 (Friedlingstein et al., 2022). Interestingly, both FREE and ASSIM capture the strengthening of the ocean carbon sink over the recent years, whereas GCP models do not. Nonetheless, it is difficult to identify an impact of the nudging on the simulated trends in ocean carbon sink in ASSIM with respect to FREE. In particular, the nudging does not improve the representation of the decadal swing of the ocean carbon sink observed before and after the 2000s. Indeed, all models' configurations fail at capturing the slowdown of the ocean carbon sink in the 2000s, including ASSIM. However, the GCP models are closer to those data-products displaying a weaker variability.

Figure 11b helps to identify the added value of the nudging by scrutinizing its impact on the simulated variability in terms of magnitude and chronology. The nudging improves the consistency between modelled and observed chronology in ocean fugacity, and slightly reinforces the magnitude of the ocean carbon sink variability. This improvement is due to the fact that fCO₂ is driven by changes in temperature and salinity in the ocean, which are directly impacted by the nudging approach. Although small, the improvement in the modelled chronology of the ocean carbon sink variability has the potential to improve the capability of the model to predict year-to-year variation in ocean carbon sink.

# 5. Skill Assessment of Key Climate and Biogeochemical Fields

## 5.1. Seasonal Timescale

ENSO diversity is considered to assess forecast performance of C3PS considering that central and eastern equatorial Pacific variability modes convey different tropical teleconnections outside the tropical Pacific (Taschetto et al., 2016; Yeh et al., 2009). For that, the forecast members are projected on the spatial patterns of the two ENSO modes shown in Figure 5 to obtain the E and C indices. As a reminder the E and C indices are uncorrelated by construction. ACC values for the start date of 1st November show very high and significant scores for all lead times (Figures 12a and 12b). C3PS performs better than persistence for lead times greater than 6 months (i.e., summer after the initialization) for the E-mode and for all lead times for the C-mode. C3PS is more skillful at predicting central Pacific ENSO variability than eastern Pacific ENSO variability, which results from
<figure><img src="image_16.png" /><figcaption>**Figure 12.** ENSO seasonal forecast skill: (a, b) anomaly correlation coefficient skill and (c, d) rms error for the ensemble-mean (30 members) as a function of lead time for the E and C indices over the period 1960–2021 for the initialization in 1st November and compared to persistence forecasts (dotted line). Red is for ASSIM as the benchmark data (i.e., potential predictability) and the blue is for HadISST1 data as the reference. Dots indicate where the correlation is significant at the 95% level based on a t test. Persistence score (or persistence forecast) refers to the autocorrelation of the referenced data set (either HadISST1 or ASSIM).</figcaption></figure>

difficulty in predicting strong El Nino events that are of E type, a common feature of seasonal prediction systems (L'Heureux et al., 2020). In general, the central equatorial Pacific is more predictable than the eastern edge, where ENSO-related phenomena involve a sharp change in convective regime and non-linear oceanic processes, resulting in a strong positive skewness of the E index (Takahashi et al., 2011). La Niña events that are of C type tend also to be more easily predictable, in particular due to the fact that they usually take place after a strong EP El Niño that have a strong discharge of heat content making the subsequent La Niña lingered over a few years (Liu et al., 2023; Sharmila et al., 2023). At last, most coupled models, like CNRM-ESM2.1, also have a warm bias in the far eastern Pacific that is influential on the forecasts (L'Heureux et al., 2022).

C3PS is effective at predicting both types of ENSO since ACC remains significant for 12 months after initialization for the E-mode, and at least 13 months for the C-mode. Potential predictability is slightly higher, the difference with ACC computed from observations increases at longer lead times. We have checked that the C3PS performances at predicting ENSO during the period 1960–2021 are comparable to the current seasonal predictions systems such as SEAS5-20C (Sharmila et al., 2023; Weisheimer et al., 2021) (not shown). RMSE scores, which take into account the prediction of the ENSO amplitude, beat persistence scores for longer lead times (Figure 12cd). Again, potential predictability is higher, in particular for ENSO-C.

ENSO skill for the 1st May starting date presents ACC values above persistence from 2 months onwards lead times (Figure 13). Again C3PS achieves better performance for the ENSO-C mode. Since boreal spring (the season of initialization) corresponds to that of the ENSO onset and the usually enhanced Madden-Julian oscillation variance, we may expect that the system has also good performance in predicting the tropical Pacific teleconnection at that season. Potential predictability exceeds that of FORCED by an average of 0.1 correlation across all lead times.
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Same as Figure 12 but for the initialization in May.</figcaption></figure>

## 5.2. Multiannual Timescale

Skill maps of ACC computed between PRED and FORCED for surface temperature show high and significant skill over large portions of the globe for forecast years Y1, Y2 and Y1-5 (Figure 14, left column). ACC scores are usually higher in the tropics than in the extra-tropics. For Y2, the skill rapidly decreases over the Eastern Pacific and the Southern Ocean but remains high and significant over the North Atlantic and Indian Oceans, Europe, Northern Asia, Northern Africa, North America and some areas in South America. Considering Y1-5, the temporal average over the five forecast years, the ACC skill is high and significant over a great portion of the Northern Hemisphere and the Indian Ocean. Potential predictability, measured by the ACC between PRED and ASSIM, is clearly higher (middle column) in many regions of the globe, including most of the continental areas, except India.

When we compare PRED and FREE skills in terms of potential predictability (Figure 14, right column), results show that some regions exhibit larger skill scores in PRED at Y1, indicating that initialization largely improves ACC scores in most of the Pacific Ocean, SPNA, western tropical Atlantic and northern South-America, central Indian Ocean and eastern Australia. In general, from Y2 onwards much of the skill is provided by the large externally-forced trend as shown by the similarity between PRED and FREE skill scores. The regions where initialization still plays an important role are the SPNA, Equatorial Pacific, Southern Pacific and Indian Oceans, as well as over North America and Brazil. For Y1-5, the added value of initialization remains over the SPNA and Southern Pacific. The fact that one of the areas of clear benefits of the ocean initialization is the SPNA is consistent with the results reported by current decadal prediction systems (IPCC, 2023).

As indicated by MSSS scores of TPI, PRED is more accurate than FREE at Y1 (Figure 15a), skill computed from both references FORCED and ASSIM showing similar scores. Comparison with observational estimates are also shown in Figure 15. After Y1, PRED and FREE performances both become statistically insignificant (Figures 14 and 15a). Similarly, focusing on the SST skill over the SPNA area, MSSS indicates that PRED performs better than FREE (Figure 15b) for lead times up to 3 years. If we focus on potential predictability, PRED is always more accurate than FREE up to Y4 over the SPNA (Figures 14 and 15c).
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Left column: anomaly correlation coefficient (ACC) skill scores for Sea Surface Temperature (SST) and surface air temperature (SAT) over land computed between PRED and FORCED for SST and between PRED and JRA55do reanalysis for SAT for forecast year 1 (top), forecast year 2 (middle) and average of forecast years (bottom). Middle column: the same but the ACC is computed between PRED and ASSIM (potential predictability). Right column: Differences between the ACC of PRED versus FREE when ASSIM is used as reference. All the data were interpolated to a regular 5-degree grid before the analysis. Stippling with gray dots indicates skill scores that are not significant at the 10% level based on block-bootstrapping as explained in the text. Stippling with light brown dots indicates ACC differences that are not significant at the 10% level based on block-bootstrapping as explained in the text.</figcaption></figure>

The impact of atmospheric initialization is also briefly addressed (see Figure S2 in Supporting Information S1). It is important to note that the atmosphere is initialized from the ERA5 reanalysis for the dynamical atmospheric variables and from ASSIM for the others. We analyzed a second set of forecasts (initialized solely on November

<figure><img src="image_19.png" /><figcaption>**Figure 15.** (a) Mean Square Skill Score (MSSS) skill scores for the SSTs for the Tripole Pacific Index (TPI). To compute MSSS, mean-square errors of PRED and FREE are compared with respect to the same reference: observations (HadISST1, blue line), FORCED (black line) and ASSIM (red line). Positive MSSS indicates that PRED performs better than FREE. The dots indicate where MSSS is statistically significant at the 10% level based on block-bootstrapping as explained in the text. (b) The same as (a) but for SSTs over the SPNA box. The TPI index is computed from raw data according to Henley et al. (2015). The SPNA index from raw data is obtained according to Bilbao et al. (2021) (SPNA: 50–65 N, 60–10 W).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Left column: anomaly correlation coefficient (ACC) skill scores for Net primary production computed between PRED and FORCED for forecast year 1 (top), forecast year 2 (middle) and average of forecast years (bottom). Middle column: the same but the ACC is computed between PRED and ASSIM (potential predictability). Right column: Differences between the ACC of PRED versus FREE when ASSIM is used as reference. All the data were interpolated to a regular 5-degree grid before the analysis. Stippling with gray dots indicates skill scores that are not significant at the 10% level based on block-bootstrapping as explained in the text. Stippling with light brown dots indicates ACC differences that are not significant at the 10% level based on block-bootstrapping as explained in the text.</figcaption></figure>

1st) using raw atmospheric initial conditions from ASSIM, without replacing the ERA5 dynamical fields in the ARPEGE restarts. Ten ensemble members were generated with this configuration. Our findings indicate that using ERA5 data for atmospheric initialization results in a slight increase in initial shock over the tropical Pacific. However, this initial shock does not significantly affect forecast skill on multi-year timescales. Additionally, atmospheric initialization moderately improves skill at seasonal timescales over the Equatorial Pacific, though we remain cautious about these results due to the relatively small ensemble size.

The forecasting skills of C3PS for biogeochemical variables such as NPP and ocean carbon fluxes are also assessed using the concept of potential predictability and compared with the skill calculated using FORCED. NPP skill scores show in general a high level of predictability over midlatitudes at Y1, Y2 and the average Y1-5 (Figure 16, left column). In contrast, the predictability of NPP in most of the tropics is very low and even the skill can be even negative when the variability of the NPP is opposite in phase with that of the target (FORCED or ASSIM). Such global features of the C3PS predictive skill for NPP contrast with the results of Séférían et al. (2014) using the IPSL-CM5A-LR model and SST anomaly initialization scheme but are in line with the findings of Frölicher et al. (2020) using GFDL-ESM2-M.

Potential predictability of NPP shows good skill scores worldwide for Y1 (Figure 16, middle column). At Y2 NPP skill decreases over some areas in the Equatorial Pacific, western Atlantic, Northern Indian and Southern Oceans, but in general it remains high and statistically significant over most of the ocean for Y1-5. Most importantly, NPP skill is high in the areas of highest marine productivity, such as the equatorial and eastern boundary upwelling systems, in particular the Canary Upwelling System. The impact of model initialization is more important on the NPP than on the SST beyond the second year of forecasting, indicating that the initialization of the BGC undoubtedly leads to benefits in predictive ability. PRED performs better than FREE practically everywhere at Y1 (Figure 16, right column). At longer horizons, ACC differences show that PRED is more accurate than FREE over
<figure><img src="image_21.png" /><figcaption>**Figure 17.** Left column: Anomaly correlation coefficient (ACC) skill scores for ocean carbon fluxes computed between PRED and FORCED for forecast year 1 (top), forecast year 2 (middle) and average of forecast years (bottom). Middle column: the same but the ACC is computed between PRED and ASSIM (potential predictability). Right column: Differences between the ACC of PRED versus FREE when ASSIM is used as reference. All the data were interpolated to a regular 5-degree grid before the analysis. Stippling with red dots indicates skill scores that are not significant at the 10% level based on block-bootstrapping as explained in the text. Stippling with gray dots indicates skill scores that are not significant at the 10% level based on block-bootstrapping as explained in the text. Stippling with light brown dots indicates ACC differences that are not significant at the 10% level based on block-bootstrapping as explained in the text.</figcaption></figure>

the Eastern North Atlantic, including the Canary Upwelling area, Tropical Atlantic, North Pacific and Central Equatorial Pacific and most of the Indian Ocean.

ACC skill of ocean carbon fluxes with FORCED as the reference is relatively high and significant over the tropical band and Southern Oceans during the first 2 years after the initialization (Figure 17, first and second column). This result is consistent with the first multi-model assessment of the ocean carbon sink prediction skills (Ilyina et al., 2021).

C3PS provides skillful predictions of ocean carbon uptake at multiannual scale over the high latitude oceans and the tropics. Potential predictability is even higher and indicates the important fact that ocean carbon fluxes can be predictable several years in advance over the areas of large carbon uptake variability such as North Atlantic and North Pacific oceans and Southern Ocean. These results support previous predictability studies based on perfect model frameworks or decadal predictions with ESMs (Lovenduski et al., 2019; Séférían et al., 2019). More importantly, this potential predictability exceeds that inferred by the knowledge of the external forcing for the first 2 years after the initialization (Figure 17, right column). After that time horizon most of the predictive skill comes from the increase of atmospheric CO2 as the primary driver of the ocean carbon sink. Within the lead years 1, 2 and 1–5, the predictable fraction of the ocean carbon sink is 37%, 19%, 16%. The fractions indicate the proportion of ocean grid points where the skill of PRED is both greater than that of FREE and statistically significant, relative to the total number of grid points in the ocean. At Y1, predictable regions include the North Atlantic and the Southern Ocean, the two major ocean carbon sink locations, as well as the Equatorial Pacific and northern Indian ocean. After Y2, only the Southern ocean carbon sink remains predictable as well as a smaller fraction of the Equatorial Pacific and Indian oceans. This result is in line with previous work made with other modeling prediction platforms (Lovenduski et al., 2019; Séférían et al., 2019).
## 6. Conclusions

In this study, the new climate prediction prototype of the CNRM-Cerfacs modeling group, C3PS is presented and evaluated. The two main novelties are that C3PS is based on an ESM, CNRM-ESM2.1, and has been designed to produce predictions from seasonal to multiannual scales. C3PS is the result of the joint work of experts in seasonal and decadal forecasting and modellers of ocean physics and biogeochemistry within the CNRM-Cerfacs research group. In addition, for multiannual predictions, C3PS has participated in the international DCPP-A exercise, and a subset of the variables produced are published in the ESGF.

The initialization procedure of C3PS consists of a full-field initialization in which all the model components are initialized from an in-house reanalysis product obtained in two steps. The first step is a forced experiment in which ocean and biogeochemistry models are driven by JRA55do reanalysis following the GCP protocol. In the second step, the T and S of this forced experiment are used to constrain only the ocean physics of CNRM-ESM2.1 through sea surface restoring and a Newtonian damping in the ocean subsurface, as described in Sanchez-Gomez et al. (2016). This method has been implemented in other climate prediction systems as in Bilbao et al. (2021). The reconstruction obtained is called dcppA-assim according to the nomenclature used in the DCPP protocol.

In this paper we have performed a basic validation of the dcppA-assim (ASSIM) experiment, which is not often done in other studies presenting forecasting systems. For us it is important to evaluate and to document the quality of our initial conditions and to investigate how the nudging of T and S affects the behavior of other variables, such as AMOC and biogeochemistry. We believe that the study of the reconstructions created to initialize the climate prediction systems is relevant, and even more so in the context of the new applications offered in the prediction of marine biogeochemistry and carbon fluxes.

ASSIM shows improvements with respect to the historical ensemble FREE in the modes of variability at the seasonal and decadal scales. The improvements are notable in the Pacific, with better representation of ENSO diversity by ASSIM and of Pacific decadal variability associated with the IPV. For other variables and other regions such as in the SPNA, ASSIM shows consistency with the temporal phase of observations in both ocean surface and heat content.

Regarding the initialization of biogeochemistry, we found the interesting result that the nudging of T and S is not sufficient to constrain the biogeochemistry, as seen in the biases presented by ASSIM in chlorophyll. We suggest that biases in nutrients, such as NO3, and an underestimation of MLD can consistently explain the misrepresentation of chlorophyll in ASSIM. This result offers perspectives for improving the reconstruction of biogeochemical variables, indicating that we should pay special attention to nutrients, which leads us to think of a nutrient nudging complementary to the nudging of physical variables.

Nevertheless we show that the T/S nudging leads to a significant improvement in the amplitude of the variability and temporal chronology of the NPP in the Tropical Pacific, coherent with previous studies (Séférrian et al., 2014). Moreover, our results also show an added value of nudging in representing carbon sink variability in terms of magnitude and timing. This improvement is due to the fact that the fugacity is controlled by changes in T and S in the ocean, which are directly affected by the nudging.

In terms of skill at seasonal scale, C3PS shows a very similar ENSO prediction skill to other seasonal forecasting systems. Considering the diversity of ENSO, the C-ENSO mode exhibits higher and significant skill levels compared to the E-ENSO mode. This is somewhat expected since the E-ENSO mode is associated with the variability of extreme El Niño events whose onsets are difficult to predict due to their nonlinear dynamics. Seasonal prediction systems also exhibit a persistent mean temperature bias in the far eastern Pacific, which alters key ENSO processes (e.g., thermocline feedback and atmospheric convection) in this region explaining the lower skill in terms of the E-ENSO mode. However the C3PS skill at seasonal timescales in the tropical Pacific is encouraging for addressing seasonal forecast skill over other regions assuming a realistic simulation of ENSO atmospheric teleconnections. Such an estimation may however suffer from the limitation of using only 30 members for the first prototype of C3PS. We have considered extending the ensemble size for future applications and evaluation.

On a multiannual scale, the C3PS results are consistent with those found in other decadal forecasting systems (IPCC, 2023), that is, C3PS shows a clear added value of ocean initialization in the prediction of SST and SAT in
the first 2 years. The novelty is a significant prediction skill of SSTs in the equatorial Pacific at Y1. On longer time scales, the added value of initialization is only detectable in the SPNA area.

The most innovative aspect of the C3PS results is the potential predictive skill displayed for NPP and carbon fluxes at different lead times. The high levels of NPP potential predictability at multiannual timescales were already reported in Séférían et al. (2014) and recently addressed in Yeager et al. (2022). These results corroborate previous findings and confirm the potential benefits for marine ecosystem prediction based on integrated physical-biogeochemical forecasting platforms such as C3PS (Tommasi et al., 2017). The fact that the evolution of carbon fluxes is potentially predictable over the regions of major carbon sink locations is also promising for improving our estimations of the future global carbon budget in the climate system.

We might consider whether running the CNRM-ESM2.1 in emission mode could affect the CO2 predictability. As demonstrated by Hajima et al. (2024, in discussion), switching from a concentration-driven to an emission-driven configuration can have direct implications for the ocean carbon sink. Testing C3PS performance in emission-driven mode is a very obvious perspective, especially since this approach is already used for natural aerosols.

To finish, although C3PS is designed with an improved initialization scheme, the C3PS multiannual predictions still suffer from initial shocks and drifts after the initialization. In particular the ENSO drift documented in Sanchez-Gomez et al. (2016) is still present in the C3PS predictions. As shown in this study, the first year after the initialization is characterized by a quasi-systematic excitation of ENSO warm events that trigger teleconnection patterns over the midlatitudes, potentially polluting the signals to be predicted. The drift problem is one of the major challenges in decadal prediction. Although some progress has been achieved since the early 2000s, drifts are still present in decadal prediction systems. In addition to improving climate models to reduce errors, another essential aspect is the improvement of the data assimilation technique to obtain initial states compatible with the climate model that will be used to make the prediction. On this line, some decadal forecasting centers are opting for “in-house reanalysis” built from the coupled models used to make the forecasts to maintain physical consistency amongst all model components. This idea involves using more complex data assimilation methods, such as the Ensemble Kalman Filter or the particle filter approaches (Counillon et al., 2014; Dai et al., 2020; Zunz et al., 2015) which have been successfully applied in the context of decadal prediction. This offers interesting pathways for improving initialization and for these reasons the implementation of a particle filter in C3PS is one of the perspectives to improve the initialization procedure.

## Data Availability Statement

The HadISST (Rayner et al., 2003) data set for SSTs and corresponding documentation are available on https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html. The EN4 ocean temperature data (Good et al., 2013) are available on https://www.metoffice.gov.uk/hadobs/en4/download.html. The RAPID array data (Moat et al., 2022) are available on https://rapid.ac.uk/challenge/data_download.php. The ESA-OC-CC data are available on https://climate.esa.int/en/projects/ocean-colour/. The simulations used in this study are performed by the CNRM-ESM2.1 ESM (Séférían, 2019) https://doi.org/10.22033/ESGF/CMIP6.1395. The model outputs (Séférían et al., 2019) are available for download on ESGF under CMIP6 projects (https://esgf-metagrid.cloud.dkrz.de/search/cmip6-dkrz/?mip_era=CMIP6&activity_id=ScenarioMIP&institution_id=CNRMCERFACS&source_id=CNRM-ESM2-1). The SURFEX-CTRP code (Decharme et al., 2019) is available at (http://www.umr-cnrm.fr/surfex/spip.php?article387) at the SURFEX website (http://www.umr-cnrm.fr/surfex). NEMO-GELATO-PISCESv2-gas (Madec et al., 2023) is also available at https://www.nemo-ocean.eu/. The access to the Git repository is granted upon request to the corresponding author. OASIS3-MCT software (Craig et al., 2017) can be downloaded at this website (https://oasis.cerfacs.fr/en/download-oasis3-mct-sources/). XIOS (Meurdesoif, 2018; Meurdesoif, 2018) can be downloaded at the XIOSV2.0 website (https://zenodo.org/records/4905653) and documentation can be found at (https://forge.ipsl.jussieu.fr/ioserver). For the ARPEGE-Climat_v6.3 code (Roehrig et al., 2020) and exact version applied to each component, please contact the authors. Most of the computations performed in this study have been done by using NCL free software version 6 (NCL 2019; https://www.ncl.ucar.edu/). Finally, a number of analyzing tools developed at CNRM, or in collaboration with CNRM scientists, is available on as Open Source code (see https://opensource.cnrm-game-meteo.fr/).
Acknowledgments

The C3PS platform has been developed in the framework of the EU-H2020 TRIATLAS project framework (https:// triatlas.w.uib.no/), under the Grant agreement 817578. The authors would like to extend their special thanks to the entire CNRM-Cerfacs team. Their contributions to the development of the modeling platform and the climate models have been instrumental in the completion of this work. We also thank the Météo-France HPC team for allowing us enough computing hours and resources to complete all the experiments in the Belenos supercomputer and Hendrix server. Boris Dewitte acknowledges support from Agencia Nacional de Investigación y Desarrollo (Concurso de Fortalecimiento al Desarrollo Científico de Centros Regionales 2020-R20F0008-CEAZA, Anillo Eclipse ACT210071, Fondecyt Regular 1231174) and the EU H2020 FutureMares project (Theme LC-CLA-06-2019, Grant agreement No 869300). Finally, the authors would like to express their sincere thanks to the anonymous reviewers for the time they spent reviewing the manuscript. Their comments and suggestions have helped us improve the quality of this study.

References

Athanasiadis, P. J., Yeager, S., Kwon, Y. O., Bellucci, A., Smith, D. W., & Tibaldi, S. (2020). Decadal predictability of North Atlantic blocking and the NAO. *npj Climate and Atmospheric Science*, 3(1), 20. https://doi.org/10.1038/s41612-020-0120

Aumont, O., Ethé, C., Tagliabue, A., Bopp, L., & Gehlen, M. (2015). PISCES-v2: An ocean biogeochemical model for carbon and ecosystem studies. *Geoscientific Model Development*, 8(8), 2465–2513. https://doi.org/10.5194/gmd-8-2465-2015

Barnier, B., Madec, G., Penduff, T., Molines, J. M., Treguier, A. M., Le Sommer, J., et al. (2006). Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution. *Ocean Dynamics*, 56(5-6), 543–567. https://doi.org/10.1007/s10236-006-0082-1

Bakker, D. C. E., Alin, S. R., Becker, M., Bittig, H. C., Castaño-Primo, R., Feely, R. A., et al. (2022). Surface ocean CO₂ atlas database version 2022 (SOCATv2022) [NCEI accession 0253659] [Dataset]. NOAA National Centers for Environmental Information. https://doi.org/10.25921/1h9f-nb73

Barnier, B., Siefridt, L., & Marchesiello, P. (1995). Thermal forcing for a global ocean circulation model using a three-year climatology of ECMWF analyses. *Journal of Marine Systems*, 6(4), 363–380. https://doi.org/10.1016/0924-7963(94)00034-9

Batté, L., & Déqué, M. (2012). A stochastic method for improving seasonal predictions. *Geophysical Research Letters*, 39(9), L09707. https://doi.org/10.1029/2012GL051406

Bellucci, A., Haarsma, R., Gualdi, S., Athanasiadis, P. J., Caian, M., Cassou, C., et al. (2015). An assessment of a multi-model ensemble of decadal climate predictions. *Climate Dynamics*, 44(9), 2787–2806. https://doi.org/10.1007/s00382-014-2164-y

Bilbao, R., Wild, S., Ortega, P., Acosta-Navarro, J., Arsouze, T., Bretonnière, P.-A., et al. (2021). Assessment of a full-field initialized decadal climate prediction system with the CMIP6 version of EC-Earth. *Earth System Dynamic*, 12(1), 173–196. https://doi.org/10.5194/esd-12-173-2021

Boer, G. J., Smith, D. M., Cassou, C., Doblas-Reyes, F., Danaba- soglu, G., Kirtman, B., et al. (2016). The decadal climate prediction project (DCPP) contribution to CMIP6. *Geoscientific Model Development*, 9(10), 3751–3777. https://doi.org/10.5194/gmd-9-3751-2016

Bojovic, D., Bilbao, R., Díaz, L. B., Donat, M., Ortega, P., Ruprich-Robert, Y., et al. (2019). The biggest unknowns related to decadal prediction: What 50 experts think are the 5 major knowledge gaps. *Bulletin of the American Meteorological Society*, 100, ES255–ES259. https://doi.org/10.1175/BAMS-D-19-0190.1

Borchert, L. F., Menary, M. B., Swingedouw, D., Sgubin, G., Hermanson, L., & Mignot, J. (2021). Improved decadal predictions of North Atlantic subpolar gyre SST in CMIP6. *Geophysical Research Letters*, 48(3), e2020GL091307. https://doi.org/10.1029/2020GL091307

Cai, W., Wang, G., Dewitte, B., Wu, L., Santoso, A., Takahashi, K., et al. (2018). Increased variability of eastern Pacific El Niño under greenhouse warming. *Nature*, 564(7735), 201–206. https://doi.org/10.1038/s41586-018-0776-9

Capotondi, A., Wittenberg, A., Kug, J. S., Takahashi, K., & McPhaden, M. (2020). ENSO diversity. In M. J. McPhaden, A. Santoso, & W. Cai (Eds.), *El Niño Southern Oscillation in a Changing Climate* (pp. 65–86). https://doi.org/10.1002/978119548164-ch4

Cassou, C., Kushnir, Y., Hawkins, E., Pirani, A., Kucharski, F., Kang, I., & Caltabiano, N. (2018). Decadal climate variability and predictability: Challenges and opportunities. *Bulletin of the American Meteorological Society*, 99(3), 479–490. https://doi.org/10.1175/BAMS-D-16-0286.1

Choi, J., & Son, S. W. (2022). Seasonal-to-decadal prediction of El Niño–Southern Oscillation and Pacific decadal Oscillation. *NPJ Climate and Atmospheric Science*, 5(1), 29. https://doi.org/10.1038/s41612-022-00251-9

Counillon, F., Bethke, I., Keenlyside, N., Bentsen, M., Bertino, L., & Zheng, F. (2014). Seasonal-to-decadal predictions with the ensemble Kalman filter and the Norwegian Earth System Model: A twin experiment. *Tellus A: Dynamic Meteorology and Oceanography*, 66(1), 21074. https://doi.org/10.3402/tellusa.v66.21074

Copernicus Climate Change Service, Climate Data Store. (2021). ORAS5 global ocean reanalysis monthly data from 1958 to present. *Copernicus Climate Change Service (C3S) Climate Data Store (CDS)*. https://doi.org/10.24381/cds.67e8eeb7

Craig, A., Valcke, S., & Coquart, L. (2017). Development and performance of a new version of the OASIS coupler, OASIS3-MCT_3.0. *Geoscientific Model Development*, 10(9), 3297–3308. https://doi.org/10.5194/gmd-10-3297-2017

Dai, P., Gao, Y., Counillon, F., Wang, Y., Kimmritz, M., & Langehaug, H. R. (2020). Seasonal to decadal predictions of regional Arctic sea ice by assimilating sea surface temperature in the Norwegian Climate Prediction Model. *Climate Dynamics*, 54(9-10), 3863–3878. https://doi.org/10.1007/s00382-020-05196-4

Delire, C., Séférian, R., Decharme, B., Alkama, R., Calvet, J.-C., Carrer, D., et al. (2020). The global land carbon cycle simulated with ISBA-CTRIP: Improvements over the last decade. *Journal of Advances in Modeling Earth Systems*, 12(9), e2019MS001886. https://doi.org/10.1029/2019MS001886

Decharme, B., Delire, C., Minvielle, M., Colin, J., Vergnes, J.-P., Alias, A., et al. (2019). Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications. *Journal of Advances in Modeling Earth Systems*, 11(5), 1207–1252. https://doi.org/10.1029/2018MS001545

Delgado-Torres, C., Donat, M. G., Gonzalez-Reviriego, N., Caron, L. P., Athanasiadis, P. J., Bretonnière, P. A., et al. (2022). Multi-model forecast quality assessment of CMIP6 decadal predictions. *Journal of Climate*, 35, 4363–4382. https://doi.org/10.1175/JCLI-D-21-0811.1

Doblas-Reyes, F., Andreu-Burillo, I., Chikamoto, Y., García-Serrano, J., Guemas, V., Kimoto, M., et al. (2013). Initialized near-term regional climate change prediction. *Nature Communications*, 4(1), 1715. https://doi.org/10.1038/ncomms2704

Doney, S. C., Mitchell, K. A., Henson, S. A., Cavan, E., DeVries, T., Gruber, N., et al. (2024). Observational and numerical modeling constraints on the global ocean biological carbon pump. *Global Biogeochemical Cycles*, 38(7), e2024GB008156. https://doi.org/10.1029/2024GB008156

Drijfhout, S., van Oldenborgh, G. J., & Cimatoribus, A. (2012). Is a decline of AMOC causing the warming hole above the North Atlantic in observed and modeled warming patterns? *Journal of Climate*, 25(24), 8373–8379. https://doi.org/10.1175/JCLI-D-12-00490.1

Dunstone, N., Lockwood, J., Solaraju-Murali, B., Reinhardt, K., Tsartsali, E. E., Athanasiadis, P. J., et al. (2022). Towards useful decadal climate services. *Bulletin of the American Meteorological Society*, 103, E1705–E1719. https://doi.org/10.1175/BAMS-D-21-0190.1

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., & Taylor, K. E. (2016). Overview of the coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. *Geoscientific Model Development*, 9(5), 1937–1958. https://doi.org/10.5194/gmd-9-1937-2016

Fan, Y., & van den Dool, H. (2008). A global monthly land surface air temperature analysis for 1948-present. *Journal of Geophysical Research*, 113(D1), D01103. https://doi.org/10.1029/2007jd008470

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Gregor, L., Hauck, J., et al. (2022). Global carbon budget 2022. *Earth System Science Data*, 14(11), 4811–4900. https://doi.org/10.5194/essd-14-4811-2022

Frölicher, T. L., Ramseyer, L., Raible, C. C., Rodgers, K. B., & Dunne, J. (2020). Potential predictability of marine ecosystem drivers. *Bio-geosciences*, 17(7), 2061–2083. https://doi.org/10.5194/bg-17-2061-2020
García-Serrano, J., Doblas-Reyes, F. J., & Coelho, C. A. S. (2012). Understanding Atlantic multi-decadal variability prediction skill. *Geophysical Research Letters*, *39*(18), L18708. https://doi.org/10.1029/2012GL053283

García-Serrano, J., Guemas, V., & Doblas-Reyes, F. J. (2015). Added-value from initialization in predictions of Atlantic multi-decadal variability. *Climate Dynamics*, *44*(9-10), 2539-2555. https://doi.org/10.1007/s00382-014-2370-7

Gidden, M. J., Riahi, K., Smith, S. J., Fujimori, S., Luderer, G., Kriegler, E., et al. (2019). Global emissions pathways under different socio-economic scenarios for use in CMIP6: A dataset of harmonized emissions trajectories through the end of the century. *Geoscientific Model Development*, *12*(4), 1443-1475. https://doi.org/10.5194/gmd-12-1443-2019

Gillett, N. P., Shiogama, H., Funke, B., Hegerl, G., Knutti, R., Matthes, K., et al. (2016). The Detection and Attribution Model Intercomparison Project (DAMIP v1.0) contribution to CMIP6. *Geoscientific Model Development*, *9*(10), 3685-3697. https://doi.org/10.5194/gmd-9-3685-2016

Good, S. A., Martin, M. J., & Rayner, N. A. (2013). EN4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates. *Journal of Geophysics Research Oceans*, *118*(12), 6704-6716. https://doi.org/10.1002/2013JC009067

Goddard, L., Kumar, A., Solomon, A., Smith, D., Boer, G., Gonzalez, P., et al. (2013). A verification framework for interannual-to-decadal predictions experiments. *Climate Dynamics*, *40*(1-2), 245-272. https://doi.org/10.1007/s00382-012-1481-2

Gouretski, V., & Reseghetti, F. (2010). On depth and temperature biases in bathythermograph data: Development of a new correction scheme based on analysis of a global ocean database. *Deep-Sea Research I*, *57*, 6-833. https://doi.org/10.1016/j.dsr.2010.03.011

Gregg, W. W., & Casey, N. W. (2004). Global and regional evaluation of the SeaWiFS chlorophyll data set. *Remote Sensing Environment*, *93*(4), 463-479. https://doi.org/10.1016/j.rse.2003.12.012

Hajima, T., Kawamiya, M., Ito, A., Tachiiri, K., Jones, C., Arora, V., et al. (2024). Consistency of global carbon budget between concentration- and emission-driven historical experiments simulated by CMIP6 Earth system models and suggestion for improved simulation of CO₂ con- centration. *EGUSphere*. https://doi.org/10.5194/egusphere-2024-188

Harris, I., Jones, P. D., Osborn, T. J., & Lister, D. H. (2014). Updated high-resolution grids of monthly climatic observations – the CRU TS3.10 Dataset. *International Journal of Climatology*, *34*(3), 623–642. https://doi.org/10.1002/joc.3711

Hauck, J., Moritz, Z., Corinne, L. Q., Gruber, N., Bakker, D. C. E., Bopp, L., et al. (2020). Consistency and challenges in the ocean carbon sink estimate for the global carbon budget. *Frontiers in Marine Science*, *7*. https://doi.org/10.3389/fmsz.2020.571720

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, *146*(730), 1999–2049. https://doi.org/10.1002/qj.3803

Henley, B. J., Gergis, J., Karoly, D. J., Power, S., Kennedy, J., & Folland, C. K. (2015). A Tripole Index for the Interdecadal Pacific Oscillation. *Climate Dynamics*, *45*(11-12), 3077-3090. https://doi.org/10.1007/s00382-015-2525-1

Huang, B., Zhu, J., Marx, L., Wu, X., Kumar, A., Hu, Z.-Z., et al. (2015). Climate drift of AMOC, North Atlantic salinity and Arctic sea ice in CFSv2 decadal predictions. *Climate Dynamics*, *44*(1-2), 559–583. https://doi.org/10.1007/s00382-014-2395-y

Huang, B., Thorne, P. W., Banzon, V. F., Boyer, T., Chepurin, G., Lawrimore, J. H., et al. (2017). Extended reconstructed sea surface temperature, version 5 (ERSSTv5): Upgrades, validations, and intercomparisons. *Journal of Climate*, *30*(20), 8179–8205. https://doi.org/10.1175/JCLI-D-16-0836.1

Ilyina, T., Li, H., Spring, A., Müller, W. A., Bopp, L., Chikamoto, M. O., et al. (2021). Predictable variations of the carbon sinks and atmospheric CO₂ growth in a multi-model framework. *Geophysical Research Letters*, *48*(6), e2020GL090695. https://doi.org/10.1029/2020GL090695

IPCC. (2023). Climate change 2023: Synthesis report. In H. Lee, & J. Romero (Eds.), *Contribution of working groups I, II and III to the sixth assessment report of the intergovernmental panel on climate change core writing team* (pp. 35–115). IPCC. https://doi.org/10.59327/IPCC/AR6-9789291691647

Smith, D., Cusack, S., Colman, A., Folland, C., Harris, G., & Murphy, J. (2007b). Improved surface temperature prediction for the coming decade from a global circulation model. *Science*, *317*(5839), 796–799. https://doi.org/10.1126/science.1139540

Karamperidou, C., Jin, F. F., & Conroy, J. L. (2017). The importance of ENSO nonlinearities in tropical pacific response to external forcing. *Climate Dynamics*, *49*(7-8), 2695-2704. https://doi.org/10.1007/s00382-016-3475-y

Keenlyside, N. S., & Ba, J. (2010). Prospects for decadal climate prediction. *Wiley Interdisciplinary Reviews: Climate Change*, *1*(5), 627–635. https://doi.org/10.1002/wcc.69

Keenlyside, N. S., Latif, M., Jungclaus, J., Kornblueh, L., & Roeckner, E. (2008). Advancing decadal-scale climate prediction in the North Atlantic sector. *Nature*, *453*(7191), 84–88. https://doi.org/10.1038/nature06921

Kirtman, B., Power, S. B., Adedoyin, J. A., Boer, G. J., Bojariu, R., Camilloni, I., et al. (2013). Near-term climate CHANGE: Projections and predictability. In T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Boschung, et al., *Climate change: The physical science basis. contribution of working group I to the fifth assessment report of the intergovernmental panel on climate change*, Cambridge University Press.

Kulk, G., Platt, T., Dingle, J., Jackson, T., Jönsson, B. F., Bouman, H. A., et al. (2020). Primary production, an index of climate change in the ocean: Satellite-based estimates over two decades. *Remote Sensing*, *12*(5), 826. https://doi.org/10.3390/rs12050826

Lee, J., Planton, Y. Y., Gleckler, P. J., Sperber, K. R., Guilyardi, E., Wittenberg, A. T., & Pallotta, G. (2021). Robust evaluation of ENSO in climate models: How many ensemble members are needed? *Geophysical Research Letters*, *48*(20), e2021GL095041. https://doi.org/10.1029/2021gl095041

L'Heureux, M. L., Levine, A. F., Newman, M., Ganter, C., Luo, J. J., Tippett, M. K., & Stockdale, T. N. (2020). ENSO prediction. In *In El Niño Southern Oscillation in a changing climate* (pp. 227–246).

L'Heureux Michelle, L., Tippett Michael, K., & Wang, W. (2022). Prediction challenges from errors in tropical Pacific Sea Surface temperature trends. *Frontiers in Climate*. https://doi.org/10.3389/fclim.2022.837483

Liu, Y., Donat, M. G., England, M. H., Alexander, L. V., Hirsch, A. L., & Delgado-Torres, C. (2023). Enhanced multi-year predictability after El Niño and La Niña events. *Nature Communications*, *14*(1), 6387. https://doi.org/10.1038/s41467-023-42113-9

Lovenduski, N. S., Yeager, S. G., Lindsay, K., & Long, M. C. (2019). Predicting near-term variability in ocean carbon uptake. *Earth System Dynamics*, *10*(1), 45–57. https://doi.org/10.5194/esd-10-45-2019

Madec, G., Bourdallé-Badie, R., Bouttier, P. A., Bricaud, C., Bruciaferri, D., Calvert, D., et al. (2017). NEMO ocean engine. https://doi.org/10.5281/ZENODO.1472492

Madec, G., & The NEMO System Team. (2023). NEMO ocean engine reference manual [Software]. Zenodo. https://doi.org/10.5281/zenodo.8167700

Materia, S., Borrelli, A., Bellucci, A., Alessandri, A., Di Pietro, P., Athanasiadis, P., et al. (2014). Impact of atmosphere and land surface initial conditions on seasonal forecasts of global surface temperature. *Journal of Climate*, *27*(24), 9253–9271. https://doi.org/10.1175/JCLI-D-14-00163.1

Mathiot, P., Jenkins, A., Harris, C., & Madec, G. (2017). Explicit representation and parametrised impacts of under ice shelf seas in the z* coordinate ocean model NEMO 3.6. *Geoscientific Model Development*, *10*(7), 2849–2874. https://doi.org/10.5194/gmd-10-2849-2017
Michou, M., Nabat, P., Saint-Martin, D., Bock, J., Decharme, B., Mallet, M., et al. (2020). Present-day and historical aerosol and ozone characteristics in CNRM CMIP6 simulations. *Journal of Advances in Modeling Earth Systems*, 12(1), e2019MS001816. https://doi.org/10.1029/2019MS001816

Meehl, G. A., Goddard, L., Murphy, J., Stouffer, R. J., Boer, G., Danabasoglu, G., et al. (2009). Decadal prediction: Can it be skillful? *Bulletin of the American Meteorological Society*, 90(10), 1467–1486. https://doi.org/10.1175/2009BAMS2778.1

Meehl, G. A., Goddard, L., Boer, G., Burgman, R., Bransta- tor, G., Cassou, C., et al. (2014). Decadal climate prediction: An update from the trenches. *Bulletin of the American Meteorological Society*, 95(2), 243–267. https://doi.org/10.1175/BAMS-D-12-00241.1

Meehl, G. A., Teng, H., Smith, D., Yeager, S., Merryfield, W., Doblas-Reyes, F., & Glanville, A. A. (2022). The effects of bias, drift, and trends in calculating anomalies for evaluating skill of seasonal-to-decadal initialized climate predictions. *Climate Dynamics*, 59(11-12), 3373–3389. https://doi.org/10.1007/s00382-022-06272-7

Meinshausen, M., Vogel, E., Nauels, A., Lorbacher, K., Meinshausen, N., Etheridge, D. M., et al. (2017). Historical greenhouse gas concentrations for climate modelling (CMIP6). *Geoscientific Model Development*, 10(5), 2057–2116. https://doi.org/10.5194/gmd-10-2057-2017

Menary, M. B., Kuhlbrodt, T., Ridley, J., Andrews, M. B., Dimdore-Miles, O. B., Deshayes, J., et al. (2018). Preindustrial control simulations with HadGEM3-GC3.1 for CMIP6. *Journal of Advances in Modeling Earth Systems*, 10(12), 3049–3075. https://doi.org/10.1029/2018MS001495

Meurdesoif, Y. (2018). *Xios fortran reference guide*. IPSL. http://forge.ipsl.jussieu.fr/ioсерven/svn/XIOS/trunk/doc/XIOS_reference_guide.pdf

Meurdesoif, Y. (2018). Xios version 2.0 [Software]. Zenodo. https://zenodo.org/records/4905653

Moat, B., Frajka-Williams, E., Smeed, D., Rayner, D., Johns, W. E., Baringer, M. O., et al. (2022). Atlantic meridional overturning circulation observed by the RAPID-MOCHA-WBTS (RAPID-meridional overturning circulation and heatflux array-western boundary time series) array at 26N from 2004 to 2020 (v2020.1) [Dataset]. *British Oceanographic Data Centre - Natural Environment Research Council*. https://doi.org/10.5285/e91b10af-6f0a-7fa7-e053-6c86abc05a09

Monerie, P.-A., Robson, J., Dong, B., & Dunstone, N. (2018). A role of the Atlantic Ocean in predicting summer surface air temperature over North East Asia? *Climate Dynamics*, 51(1), 473–491. https://doi.org/10.1007/s00382-017-3935-z

Muller, R. A., Curry, J., Groom, D., Jacobsen, R., Perlmutter, S., Rohde, R., et al. (2013). Decadal variations in the global atmospheric land temperatures. *Journal of Geophysical Research: Atmospheres*, 118(11), 5280–5286. https://doi.org/10.1002/jgrd.50458

O'Kane, T. J., Scaife, A. A., Kushnir, Y., Brookshaw, A., Buontempo, C., Carlin, D., et al. (2023). Recent applications and potential of near-term (interannual to decadal) climate predictions. *Frontiers in Climate*, 5. https://doi.org/10.3389/fclim.2023.1121626

O'Neill, B. C., Tebaldi, C., van Vuuren, D. P., Eyring, V., Friedlingstein, P., Hurtt, G., et al. (2016). The scenario model intercomparison project (ScenarioMIP) for CMIP6. *Geoscientific Model Development*, 9(9), 3461–3482. https://doi.org/10.5194/gmd-9-3461-2016

Park, J.-Y., Stock, C. A., Dunne, J. P., Yang, X., & Rosati, A. (2019). Seasonal to multiannual marine ecosystem prediction with a global Earth system model. *Science*, 365(6450), 284–288. https://doi.org/10.1126/science.aav6634

Pohlmann, H., Jungclaus, J. H., Köhl, A., Stammer, D., & Marotzke, J. (2009). Initializing Decadal Climate Predictions with the GECCO Oceanic Synthesis: Effects on the North Atlantic. *Journal of Climate*, 22(14), 3926–3938. https://doi.org/10.1175/2009JCLI2535.1

Pohlmann, H., Kröger, J., Greatbatch, R. J., & Müller, W. A. (2017). Initialization shock in decadal hindcasts due to errors in wind stress over the tropical Pacific. *Climate Dynamics*, 49(7-8), 2685–2693. https://doi.org/10.1007/s00382-016-3486-8

Rayner, N. A., Parker, D. E., Horton, E. B., Folland, C. K., Alexander, L. V., Rowell, D. P., et al. (2003). Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. *Journal of Geophysical Research*, 108(D14), 4407. https://doi.org/10.1029/2002JD002670

Rohde, R., Muller, R. A., Jacobsen, R., Muller, E., Perlmutter, S., et al. (2013). A new estimate of the average earth surface land temperature spanning 1753 to 2011. *Geoinfor Geostat: An Overview*, 1(1). https://doi.org/10.4172/2327-4581.1000101

Roehrig, R., Beau, I., Saint-Martin, D., Alias, A., Decharme, B., Guérémy, J.-F., et al. (2020). The CNRM global atmosphere model ARPEGE-climate 6.3: Description and evaluation. *Journal of Advances in Modeling Earth Systems*, 12(7), e2020MS002075. https://doi.org/10.1029/2020MS002075

Sanchez-Gomez, E., Cassou, C., Ruprich-Robert, Y., Fernandez, E., & Terray, L. (2016). Drift dynamics in a coupled model initialized for decadal forecasts. *Climate Dynamics*, 46(5-6), 1819–1840. https://doi.org/10.1007/s00382-015-2678-y

Salas Mélia, D. (2002). A global coupled sea ice-ocean model. *Ocean Modelling*, 4(2), 137–172. https://doi.org/10.1016/S1463-5003(01)00015-4

Schuster, M., Grieger, J., Richling, A., Schartner, T., Illing, S., Kadow, C., et al. (2019). Improvement in the decadal prediction skill of the North Atlantic extratropical winter circulation through increased model resolution. *Earth System Dynamics*, 10(4), 901–917. https://doi.org/10.5194/esd-10-901-2019

Séférian, R., Nabat, P., Michou, M., Saint-Martin, D., Voldoire, A., & Colin, J. (2019). Evaluation of CNRM earth-system model, CNRM-ESM2-1: Role of Earth system processes in present-day and future climate. *Journal of Advances in Modeling Earth Systems*, 11(12), 4182–4227. https://doi.org/10.1029/2019MS001791

Seferian, R. (2019). CNRM-CERFACS CNRM-ESM2-1 model output prepared for CMIP6 ScenarioMIP. Version YYYYMMDD [Software]. *Earth System Grid Federation*. https://doi.org/10.22033/ESGF/CMIP6.1395

Séférian, R., Berthet, S., Yool, A., Palmiéri, J., Bopp, L., Tagliabue, A., et al. (2020). Tracking improvement in simulated marine biogeochemistry between CMIP5 and CMIP6. *Current Climate Change Reports*, 6(3), 95–119. https://doi.org/10.1007/s40641-020-00160-0

Servonnat, J., Mignot, J., Guilyardi, E., Swingedouw, D., Séférian, R., & Labetoulle, S. (2015). Reconstructing the sub-surface ocean decadal variability using surface nudging in a perfect model framework. *Climate Dynamics*, 44(1-2), 315–338. https://doi.org/10.1007/s00382-014-2184-7

Séférian, R., Bopp, L., Gehlen, M., Swingedouw, D., Mignot, J., Guilyardi, E., & Servonnat, J. (2014). Multiyear predictability of tropical marine productivity. *Proceedings of the National Academy of Sciences*, 111(32), 11646–11651. https://doi.org/10.1073/pnas.1315855111

Séférian, R., Berthet, S., & Chevallier, M. (2018). Assessing the decadal predictability of land and ocean carbon uptake. *Geophysical Research Letters*, 45(5), 2455–2466. https://doi.org/10.1002/2017GL076092

Sharmila, S., Hendon, H., Alves, O., Weisheimer, A., & Balmaseda, M. (2023). Contrasting El Niño-La Niña predictability and prediction skill in 2-year reforecasts of the twentieth century. *Journal of Climate*, 36(5), 1269–1285. https://doi.org/10.1175/JCLI-D-22-0028.1

Smith, D., Cusack, S., Colman, A., Folland, C., Harris, G., & Murphy, J. (2007). Improved surface temperature prediction for the coming decade from a global circulation model. *Science*, 317(5839), 796–799. https://doi.org/10.1126/science.1139540

Smith, D. M., Eade, R., Scaife, A. A., Caron, L. P., Danabasoglu, G., DelSole, T. M., et al. (2019). Robust skill of decadal climate predictions. *NPJ Climate and Atmospheric Science*, 2(1), 13. https://doi.org/10.1038/s41612-019-0071-y

Smith, D. M., Scaife, A. A., Eade, R., Athanasiadis, P., Bellucci, A., Bethke, I., et al. (2020). North Atlantic climate far more predictable than models imply. *Nature*, 583(7818), 796–800. https://doi.org/10.1038/s41586-020-2525-0

Swingedouw, D., Bily, A., Esquerdo, C., Borchert, L. F., Sgubin, G., Mignot, J., & Menary, M. (2021). On the risk of abrupt changes in the North Atlantic subpolar gyre in CMIP6 models. *Annals of the New York Academy of Sciences*, 15(4), 187–201. https://doi.org/10.1111/nyas.14659
Takahashi, K., Montecinos, A., Goubanova, K., & Dewitte, B. (2011). ENSO regimes: Reinterpreting the canonical and Modoki El Niño. *Geophysical Research Letters*, 38(10). https://doi.org/10.1029/2011GL047364

Taschetto, A. S., Rodrigues, R. R., Meehl, G. A., McGregor, S., & England, M. H. (2016). How sensitive are the Pacific–tropical North Atlantic teleconnections to the position and intensity of El Niño-related warming? *Climate Dynamics*, 46(5-6), 1841–1860. https://doi.org/10.1007/s00382-015-2679-x

Terhaar, J., Goris, N., Müller, J. D., DeVries, T., Gruber, N., Hauck, J., et al. (2024). Assessment of global ocean biogeochemistry models for ocean carbon sink estimates in RECCAP2 and recommendations for future studies. *Journal of Advances in Modeling Earth Systems*, 16(3), e2023MS003840. https://doi.org/10.1029/2023MS003840

Tittensor, D. P., Beger, M., Boerder, K., Boyce, D. G., Cavanagh, R. D., Cosandey-Godin, A., et al. (2019). Integrating climate adaptation and biodiversity conservation in the global ocean. *Science Advances*, 5(11), eaay9969. https://doi.org/10.1126/sciadv.aay9969

Tietsche, S., Balmaseda, M., Zuo, H., Roberts, C., Mayer, M., & Ferranti, L. (2020). The importance of North Atlantic Ocean transports for seasonal forecasts. *Climate Dynamics*, 55(7-8), 1995–2011. https://doi.org/10.1007/s00382-020-05364-6

Tommasi, D., Stock, C. A., Hobday, A. J., Methot, R., Kaplan, I. C., Eveson, J. P., et al. (2017). Managing living marine resources in a dynamic environment: The role of seasonal to decadal climate forecasts. *Progress in Oceanography*, 152, 15–49. https://doi.org/10.1016/j.pocean.2016.12.011

Tsujino, H., Urakawa, L. S., Griffies, S. M., Danabasoglu, G., Adcroft, A. J., Amaral, A. E., et al. (2020). Evaluation of global ocean–sea-ice model simulations based on the experimental protocols of the Ocean Model Intercomparison Project phase 2 (OMIP-2). *Geoscientific Model Development*, 13(8), 3643–3708. https://doi.org/10.5194/gmd-13-3643-2020

Timmermann, A., An, S. I., Kug, J. S., Jin, F. F., Cai, W., Capotondi, A., et al. (2018). El Niño–Southern Oscillation complexity. *Nature*, 559(7715), 535–545. https://doi.org/10.1038/s41586-018-0252-6

The NCAR Command Language. (2019). The NCAR Command Language (Version 6.6.2) [Software]. UCAR/NCAR/CISL/TDD. https://doi.org/10.5065/D6WD3XH5

Valente, A., Sathyendranath, S., Brotas, V., Groom, S., Grant, M., Jackson, T., et al. (2022). A compilation of global bio-optical in situ data for ocean colour satellite applications – version three. *Earth System Science Data*, 14(12), 5737–5770. https://doi.org/10.5194/essd-14-5737-2022

Verfaillie, D., Doblas-Reyes, F. J., Donat, M. G., Pérez-Zanón, N., Solaraju-Murali, B., et al. (2020). How reliable are decadal climate predictions of near-surface air temperature? *Journal of Climate*, 34(2), 697–713. https://doi.org/10.1175/JCLI-D-20-0138.1

Voldoire, A., Saint-Martin, D., Sénési, S., Decharme, B., Alias, A., Chevallier, M., et al. (2019). Evaluation of CMIP6 DECK experiments with CNRM-CM6-1. *Journal of Advances in Modeling Earth Systems*, 11(7), 2177–2213. https://doi.org/10.1029/2019MS001683

Waldman, R., Hirschi, J., Voldoire, A., Cassou, C., & Msadek, R. (2021). Clarifying the Relation between AMOC and Thermal Wind: Application to the Centennial Variability in a Coupled Climate Model. *J. Phys. Oceanogr.*, 51(2), 343–364. https://doi.org/10.1175/JPO-D-19-0284.1

Wanninkhof, R. (2014). Relationship between wind speed and gas exchange over the ocean revisited. *Limnology and Oceanography: Methods*, 12(6), 351–362. https://doi.org/10.4319/lom.2014.12.351

Weisheimer, A., Balmaseda, M., Stockdale, T., Mayer, M., de Boisseson, E., Senan, R., & Johnson, S. (2021). *Retrospective two-year ENSO predictions during the 20th century* (Vol. 169, pp. 7–8). ECMWF Newsletter. https://www.ecmwfint/sites/default/files/elibrary/2021/20225-newsletter-no-169-autumn-2021.pdf

Wu, B., Zhou, T., Li, C., Müller, W. A., & Lin, J. (2019). Improved decadal prediction of Northern-Hemisphere summer land temperature. *Climate Dynamics*, 53(3), 1357–1369. https://doi.org/10.1007/s00382-019-04658-8

Yeager, S. G., Danabasoglu, G., Rosenbloom, N. A., Strand, W., Bates, S. C., Meehl, G. A., et al. (2018). Predicting near-term changes in the Earth System: A large ensemble of initialized decadal prediction simulations using the Community Earth System Model. *Bulletin of the American Meteorological Society*, 99(9), 1867–1886. https://doi.org/10.1175/bams-d-17-0098.1

Yeager, S. G., Rosenbloom, N., Glanville, A. A., Wu, X., Simpson, I., Li, H., et al. (2022). The Seasonal-to-Multiyear Large ensemble (SMYLE) prediction system using the Community Earth System Model version 2. *Geoscientific Model Development*, 15(16), 6451–6493. https://doi.org/10.5194/gmd-15-6451-2022

Yeh, S. W., Kug, J. S., Dewitte, B., Kwon, M. H., Kirtman, B. P., & Jin, F. F. (2009). El Niño in a changing climate. *Nature*, 461(7263), 511–514. https://doi.org/10.1038/nature08316

Zunz, V., Goosse, H., & Dubinkina, S. (2015). Impact of the initialization on the predictability of Southern Ocean sea ice at interannual to multi-decadal timescales. *Climate Dynamics*, 44(7-8), 2267–2286. https://doi.org/10.1007/s00382-014-2344-9

## References From the Supporting Information

Boyer, T. P., Garcia, H. E., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Reagan, J. R., et al. (2018). World Ocean Atlas 2018 [Dataset]. NOAA National Centers for Environmental Information. https://www.ncei.noaa.gov/archive/accession/NCEI-WOA18