Presentation and Evaluation of the IPSL‐CM6A‐LR Climate Model

Abstract

This study presents the global climate model IPSL‐CM6A‐LR developed at Institut Pierre‐Simon Laplace (IPSL) to study natural climate variability and climate response to natural and anthropogenic forcings as part of the sixth phase of the Coupled Model Intercomparison Project (CMIP6). This article describes the different model components, their coupling, and the simulated climate in comparison to previous model versions. We focus here on the representation of the physical climate along with the main characteristics of the global carbon cycle. The model's climatology, as assessed from a range of metrics (related in particular to radiation, temperature, precipitation, and wind), is strongly improved in comparison to previous model versions. Although they are reduced, a number of known biases and shortcomings (e.g., double Intertropical Convergence Zone [ITCZ], frequency of midlatitude wintertime blockings, and El Niño–Southern Oscillation [ENSO] dynamics) persist. The equilibrium climate sensitivity and transient climate response have both increased from the previous climate model IPSL‐CM5A‐LR used in CMIP5. A large ensemble of more than 30 members for the historical period (1850–2018) and a smaller ensemble for a range of emissions scenarios (until 2100 and 2300) are also presented and discussed.

Introduction

The Institut Pierre-Simon Laplace Climate Modelling Centre (IPSL CMC, see https://cmc.ipsl.fr) has set up a new version of its climate model in the runup of Phase 6 of the Coupled Model Intercomparison Project (known as CMIP6; see Eyring et al., 2016, for more information). Here we provide a brief description of the coupled model, document the model climatology and its performance against a range of observations and reanalyses, and present some key emerging properties of the model (internal variability and response to forcings). The implementation of the model boundary conditions (Lurton et al., 2020) and the development process for this new model configuration in preparation to CMIP6 are described in two companion papers.

IPSL CMC developed IPSL-CM5A-LR (CM stands for climate model and LR for low resolution) as its main model for Phase 5 of CMIP (Dufresne et al., 2013; Szopa et al., 2013). IPSL-CM5A-LR also had two variants: a medium-resolution configuration, IPSL-CM5A-MR, and an experimental version, IPSL-CM5B-LR, based on a new version of the atmospheric physics (Hourdin et al., 2013a). The resolution of the atmospheric model was 96 × 95 points in longitude and latitude in the LR configuration, and 144 × 143 in the MR configuration. Both versions had 39 layers in the vertical. The nominal resolution of the NEMO oceanic model was 2° for both configurations. Since then, many improvements have been implemented in the various model components: LMDZ (atmosphere), NEMO (ocean, sea ice, marine biogeochemistry), and ORCHIDEE (land surface, hydrology, land carbon cycle). In this article we describe only the coupled ocean-atmosphere model and the carbon cycle in the terrestrial and marine model components as the full Earth System version of IPSL-CM6 is still under development. The resolution of the atmospheric model is now 144 × 143 points in longitude and latitude, which corresponds to an average resolution of $\sqrt{(}4\pi R^2/144/142)=157$ km (R being the Earth's radius) and 79 vertical layers (with a model top at ∼80 km). The low horizontal resolution (LR) of IPSL-CM6 thus corresponds to the medium horizontal resolution (MR) of IPSL-CM5. The nominal resolution of the ocean model has been increased to 1° and 75 layers in the vertical.

This article provides an entry point to the IPSL-CM6A-LR model with a brief scientific and technical description of the model, a thorough evaluation of its climatology, and some presentation of the DECK (Diagnostic, Evaluation, and Characterization of Klima) and ScenarioMIP simulations that were prepared for CMIP6. Further studies on the IPSL-CM6A-LR model emerging properties and model intercomparison studies are expected to be ongoing for the next few years.

Brief Overview of the IPSL-CM6A-LR Model Introduction

IPSL-CM6A-LR is composed of the LMDZ atmospheric model Version 6A-LR (Hourdin et al., 2020a), the NEMO oceanic model Version 3.6 (see references below) and the ORCHIDEE land surface model Version 2.0. We briefly describe below each of the three model components and the coupling procedure between them. Further description of the IPSL-CM6A-LR climate model is available on the ES-DOC interface (https://explore.es-doc.org/cmip6/models/ipsl/ipsl-cm6a-lr).

LMDZ6 Atmospheric Component

The atmospheric general circulation model LMDZ6A-LR is based on a finite-difference formulation of the primitive equations of meteorology (Sadourny & Laval, 1984), on a staggered and stretchable longitude-latitude grid (the Z of LMDZ standing for Zoom). Water vapor, liquid and solid water, and atmospheric trace species are advected with a monotonic second-order finite volume scheme (Hourdin & Armengaud, 1999; Van Leer, 1977). In the vertical, the model uses a classic so-called hybrid sigma-pressure coordinate. Regarding the physical parameterizations, IPSL participated in CMIP5 with two versions: a “Standard Physics” version (atmospheric component LMDZ5A used in IPSL-CM5A Hourdin et al., 2013b) and a “New Physics” (NP) version (LMDZ5B used in IPSL-CM5B Hourdin et al., 2013b) based on a full rethinking of the parameterizations of turbulence, convection, and clouds on which the 6A version is built. This NP package includes in particular a turbulent scheme based on the prognostic equation for the turbulent kinetic energy that follows Yamada (1983), a mass flux representation of the organized structures of the convective boundary layer called “Thermal Plume Model” (Hourdin et al., 2002; Rio & Hourdin, 2008; Rio et al., 2010), and a parameterization of the cold pools or wakes created below cumulonimbus by the evaporation of convective rainfall (Grandpeix & Lafore, 2010; Grandpeix et al., 2010). The “episodic mixing and buoyancy sorting” scheme originally developed by Emanuel (1991) used for deep convection was modified to make the closure and triggering rely on the description of the subcloud vertical motions by thermal plumes and wakes (Rio et al., 2009). Regarding convection, two important improvements were made from Versions 5B to 6A: a modification of the lateral detrainment in the thermal plume model that allows to represent satisfactorily well the transition from stratocumulus to cumulus clouds (Hourdin et al., 2019a) and the introduction of a statistical triggering for deep convection (Rochetin et al., 2014a, 2014b). The radiation scheme was inherited from the European Centre for Medium-Range Weather Forecasts. In the LMDZ6A version, it includes the Rapid Radiative Transfer Model (RRTM) code for thermal infrared radiation and an improved six-band version of the Fouquart and Bonnel (2014b) scheme for solar radiation. Cloud cover and cloud water content are computed using a statistical scheme using a lognormal function for deep convection (Bony & Emanuel, 2001) and a bigaussian function for shallow cumulus (Jam et al., 2013).

The 6A-LR version is based on a regular horizontal grid with 144 points regularly spaced in longitude and 142 in latitude, corresponding to a resolution of 2.5° × 1.3°. The model has 79 vertical layers and extends up to 80 km, which makes it a “high-top” model. It includes a representation of gravity waves generated by mountains as well as by convection (Lott & Guez, 2013) and fronts (de la Cámara & Lott, 2015; de la Cámara et al., 2016). The model shows a self-generated quasi-biennal oscillation (QBO) whose period has been tuned to the observed one for the present-day climate. The source of water vapor in the stratosphere due to methane oxidation is not activated.

The readers are directed to Hourdin, Rio, Grandpeix, et al. (2020a) for more details.

NEMO Oceanic Component

The ocean component used for IPSL-CM6A-LR is based on the Version 3.6 stable of NEMO (Nucleus for European Models of the Ocean), which includes three major components: the ocean physics NEMO-OPA (Madec et al., 2017), the sea ice dynamics and thermodynamics NEMO-LIM3 (Rousset et al., 2015; Vancoppenolle et al., 2009), and the ocean biogeochemistry NEMO-PISCES (Aumont et al., 2015). The configuration used is eORCA1 (with the e standing for extended), the quasi-isotropic global tripolar grid with a 1° nominal resolution, and extended to the south so as to better represent the contribution of Antarctic under-ice shelf seas to the Southern Ocean freshwater cycle (Mathiot et al., 2017). The grid has a latitudinal grid refinement of 1/3° in the equatorial region. Vertical discretization uses a partial step formulation (Barnier et al., 2006), which ensures a better representation of bottom bathymetry, with 75 levels. The initial layer thicknesses increase nonuniformly from 1 m at the surface to 10 m at 100 m depth and reaches 200 m at the bottom; they are subsequently time dependent (Levier et al., 2007).

Ocean Physics: NEMO-OPA

The eORCA1 configuration used has a nonlinear free surface using the variable volume layer formulation, which induces time variability of all layer thicknesses (Levier et al., 2007). It uses a polynomial representation of the equation of state TEOS-10 (Roquet et al., 2015). The vertical mixing of tracers and momentum uses the turbulent kinetic energy scheme (Blanke & Delecluse, 1993; Gaspar et al., 1990) and an energy-constrained parameterization of mixing due to internal tides (de Lavergne, 2016; de Lavergne et al., 2019). There is no constant background diffusivity other than a floor at molecular levels: 1.4 × 10⁻⁶ m² s⁻¹ for momentum and 1.4 × 10⁻⁷ m² s⁻¹ for tracers. The mixing induced by submesoscale processes in the mixed layer is also parameterized (Fox-Kemper et al., 2011). A quadratic bottom friction boundary condition is applied together with a parameterization of a diffusive bottom boundary layer for the tracers with a coefficient of 1,000 m² s⁻¹. The model uses an energy-enstrophy-conserving scheme for momentum advection, and a no-slip boundary condition is applied on the momentum equations. Lateral diffusion of momentum is performed on geopotential surfaces and uses a Laplacian viscosity with a coefficient of 20,000 m² s⁻¹. Lateral diffusion of tracers is performed along isoneutral surfaces using Laplacian mixing with a spatially varying coefficient of 1,000 m² s⁻¹ at the equator decreasing with the reduction of the grid spacing with the latitude and reaches a value less than 500 m² s⁻¹ poleward to 60°N and 60°S. In addition, there is a parameterization of adiabatic eddy mixing (Gent & Mcwilliams, 1990) varying spatially as a function of Rossby radius and local growth rate of baroclinic instabilities. The configuration also includes representation of the interaction between incoming shortwave radiation into the ocean and the phytoplankton (Lengaigne et al., 2009). A spatially varying geothermal heat flux is applied at the bottom of the ocean (Goutorbe et al., 2011), with a global mean value of 66 mW m⁻².

Sea Ice: NEMO-LIM

IPSL-CM6A-LR utilizes v3.6 of the Louvain-la-Neuve Ice Model (LIM), instead of Version 2 for IPSL-CM5; hence, many of the sea ice model features were revised since CMIP5. LIM3.6 is a multicategory halothermodynamic dynamic sea ice model embedded in the NEMO environment (Rousset et al., 2015; Vancoppenolle et al., 2009), based on the Arctic Ice Dynamics Joint EXperiment (AIDJEX) framework (Coon et al., 1974). LIM3.6 combines the ice thickness distribution approach (Bitz et al., 2001; Lipscomb, 2001; Thorndike et al.,1975), the conservation of horizontal momentum (Hibler, 1979), treating sea ice as a 2-D elastic-viscous plastic continuum (Bouillon et al., 2013; Hunke & Dukowicz, 1997), horizontal transport (Prather, 1986), and energy-conserving halothermodynamics (Bitz & Lipscomb, 1999; Vancoppenolle et al., 2009). The multiple ice categories allow resolving the enhanced growth of thin ice and solar radiation uptake through thin ice and the redistribution of thin onto thick ice through ridging and rafting. Sea ice salinity is an integral part of the model, evolving dynamically to resolve brine entrapment and drainage and influencing sea ice thermal properties and ice-ocean exchanges (Vancoppenolle et al., 2009). Five thickness categories are used. Ice temperature and salinity fields are further discretized onto two vertical layers of sea ice and one layer of snow. Horizontally, the ice fields are resolved on the same grid as the ocean component.

The large-scale sea ice state was first adjusted in ocean-sea-ice-only simulations, and then at the end of each tuning cycle of a long fully-coupled simulation, which adds up to several thousands of years. The sea ice tuning parameters include the cloud-sky albedo nodal values (dry snow 0.87, wet snow 0.82, dry ice 0.65, and wet ice 0.58) and the snow thermal conductivity (0.31 W m⁻¹ K⁻¹). The albedo values lie in the high end of the range, to enhance sea ice formation and reduce melting, and compensate for the effects of high air temperatures above sea ice, in particular in the Arctic winter. The albedo nodal values were kept within observational uncertainty range, leaving a low Arctic sea ice bias still. The ice strength parameter was set to P^∗=20,000 N m⁻². A maximum ice concentration is imposed, which is equivalent to imposing a minimum open water fraction and done specifically for each hemisphere (0.997 in the Northern Hemisphere and 0.95 in the Southern Hemisphere). This choice is justified by the difficulty of the model to maintain open water within the pack, in particular in winter, and even more so for Antarctic sea ice.

Ocean Biogeochemistry: NEMO-PISCES

The biogeochemical model is based on PISCES-v2 (Aumont et al., 2015) which simulates the lower trophic levels of marine ecosystem (phytoplankton, microzooplankton, and mesozooplankton) and the biogeochemical cycles of carbon and of the main nutrients (P, N, Fe, and Si). There are 24 prognostic variables (tracers) including two phytoplankton compartments (diatoms and nanophytoplankton), two zooplankton size classes (microzooplankton and mesozooplankton), and a description of the carbonate chemistry. Formulations in PISCES-v2 are based on a mixed Monod/Quota formalism. On the one hand, stoichiometry of C/N/P is fixed and growth rate of phytoplankton is limited by the external availability in N, P, and Si. On the other hand, the iron and silicon quotas are variable and growth rate of phytoplankton is limited by the internal availability in Fe. Nutrients and/or carbon is supplied to the ocean from three different sources: atmospheric deposition, rivers, and sediment mobilization. PISCES is used here to compute air-sea fluxes of carbon and also the effect of a biophysical coupling: the chlorophyll concentration produced by the biological component feedbacks on the ocean heat budget by modulating the absorption of light as well as the oceanic heating rate (Lengaigne et al., 2009).

ORCHIDEE Land Surface Component

ORCHIDEE is a global process-based model of the land surface and the terrestrial biosphere that calculates water, energy, and carbon fluxes between the surface and the atmosphere. The model, initially described in Krinner et al. (2005) for the version used in the IPSL-CM5 model, has been significantly improved in Version 2.0 used in IPSL-CM6A-LR. We only summarize below the main characteristics of ORCHIDEE and key improvements from the CMIP5 version.

The vegetation heterogeneity is described using fractions of 15 different plant functional types (PFTs Prentice et al., 1992) for each grid cell. All PFTs share the same equations but with different parameters, except for the leaf phenology. The annual evolution of the PFT maps (including a wood harvest product) is derived from the LUHv2 database (Lurton et al., 2020). In each grid cell, the PFTs are grouped into three soil tiles according to their physiological behavior: high vegetation (forests) with eight PFTs, low vegetation (grasses and crops) with six PFTs, and bare soil with one PFT. An independent hydrological budget is calculated for each soil tile, to prevent forests from exhausting all soil moisture. In contrast, only one energy budget (and snow budget) is calculated for the whole grid cell. Note that the energy budget is solved with an implicit numerical scheme that couples the lower atmosphere to the surface, in order to increase numerical stability. All components of the surface energy and water budgets, as well as plant/soil carbon fluxes, are computed at the same time step as the atmospheric physics (i.e., 15 min; Hourdin, Rio, Grandpeix, et al., 2020a) using a standard “big leaf” approach, but the “slow” processes (carbon allocation in the different plant reservoir and litter and soil carbon dynamic) are computed on a daily time step. The routing scheme to transform runoff into river discharge to the ocean (Ngo-Duc et al., 2007) also proceeds at the daily time step and has not changed since IPSL-CM5.

A physically based 11-layer soil hydrology scheme has replaced the two-layer bucket model used in IPSL-CM5. Vertical water fluxes are described using the Richard equation discretized with 11 layers for a 2 m soil depth, and a free drainage condition is imposed at the bottom of the soil column (de Rosnay et al., 2002; D'Orgeval et al., 2008). As detailed in Wang et al. (2016), the vertical discretization for heat diffusion is now identical to that adopted for water up to 2 m. Furthermore, the soil depth for heat diffusion is extended to 90 m, with a zero flux condition at the bottom and 18 calculation nodes, extrapolating the moisture content of the deepest hydrological layer to the entire profile between 2 and 90 m. The soil thermal properties (heat capacity and conductivity) of each layer now depend on soil moisture and soil texture, like the soil hydrological properties (hydraulic conductivity and diffusivity). Each model grid cell is characterized by the dominant soil texture, as derived from the map of Zobler (1986) (but reduced to three classes: coarse/sandy loam, medium/loam, and fine/clay loam), and controlling the constant soil parameters (porosity, Van Genuchten parameters, field capacity and wilting point, and dry and saturated thermal properties). All these changes have a significant impact on the surface temperature and its high-frequency variability in most regions (Cheruy et al., 2017).

In contrast to IPSL-CM5, soil freezing is allowed and diagnosed in each soil layer following a scheme proposed by Gouttevin et al. (2012), but the latent heat release/consumption associated with water freezing/thawing is not accounted for. The freezing state of the soil mainly impacts the computation of soil thermal and hydraulic properties, reducing, for instance, the water infiltration capacity at soil surface. Finally, the one-layer snow scheme of IPSL-CM5 was replaced by a three-layer scheme of intermediate complexity described in Wang et al. (2013) and inspired by the scheme proposed in Boone and Etchevers (2001). A revised parameterization of the vegetation and snow albedo has been also introduced with optimized parameters based on remote sensing albedo data from MODIS sensor.

For the carbon cycle, photosynthesis depends on light availability, CO₂ concentration, soil moisture, and surface air temperature. It is parameterized based on Farquhar et al. (1980) and Collatz et al. (1992) for C3 and C4 plants, respectively. We used the implementation proposed by Yin and Struik (2009) that derives an analytical solution of the three equations linking the net assimilation rate, the stomatal conductance, and the intercellular CO₂ partial pressure. In addition, the new version of ORCHIDEE used in IPSL-CM6A-LR includes a “downregulation” capability which accounts for a reduction of the maximum photosynthesis rates as the CO₂ concentration increases in order to account for nutrient limitations. This downregulation mechanism is modeled as a logarithmic function of the CO₂ concentration relative to 380 ppm following Sellers et al. (1996). Once the carbon is fixed by photosynthesis, we compute the autotrophic respiration (growth and maintenance) and then allocate the remaining carbon into eight plant compartments (below and above ground sapwood and heartwood, leaves, fruit, roots, and reserves). Each compartment has a specific turnover depending on environmental stresses, and the living biomass is turned into a litter pool that is distributed in four compartments (metabolic or structural, both above or below ground). The litter is decomposed following first-order kinetics equations, modulated by upper soil moisture and temperature, with a fraction that is respired and a fraction that is distributed into three soil organic carbon pools (active, slow, and passive), following the CENTURY model (Parton et al., 1987). Each soil organic carbon pool is also decomposed following first-order kinetic equations modulated by soil moisture and temperature. Overall, the carbon respired from the litter and soil carbon pools defines the heterotrophic respiration.

Coupling Between the Components

The LMDZ and ORCHIDEE models are coupled at every time step of the physics of the atmospheric model (i.e., 15 min) with the exception of the biogeochemical processes and the vegetation dynamics for which the coupling frequency is 1 day.

The coupling between LMDZ and NEMO in IPSL-CM5 is described in Marti et al. (2010). It is now performed with the OASIS3-MCT coupler. IPSL-CM6A-LR introduces some modifications and new features: models are coupled with a frequency of 90 min, which is both the time step of the sea ice model and of the radiation computation in the atmosphere. Atmospheric variables passed to the ocean model (heat, water, and momentum fluxes) are averaged temporally over the six 15-min time steps of the LMDZ physics. The flux of freshwater from rivers is passed to the ocean model at river mouth locations with a frequency of 1 day, which is also the time step of the river routing in ORCHIDEE. On the ocean grid, the water coming from a river is smoothed over about 200 km to avoid strong haloclines that may occasionally cause the ocean model to crash. To ensure water conservation, the water flux into endorheic basins is globally integrated and homogeneously redistributed over the ocean. Oceanic model variables sent to the atmosphere every 90 min are sea surface temperature (SST), sea ice fraction, sea ice surface temperature, and albedo, averaged on two ocean dynamics time steps. Albedo for the open surface ocean is computed at every time step in LMDZ following Séférian et al. (2018), which represents a significant improvement over the parameterization used in IPSL-CM5 models. Ocean albedo is a function of solar zenith angle, waveband, and surface wind speed; the optional dependence on chlorophyll content of the surface ocean has not been activated. Separate albedoes are provided to the radiative transfer scheme for direct and diffuse radiation.

The model includes a very simple scheme to represent the water budget of ice sheets. Snow can accumulate on the land ice fraction of a grid box, while water vapor can deposit or sublimate depending on the surface relative humidity. The snowpack is capped to a value of 3,000 kg m⁻², and any excess is sent to a buffer reservoir before returning to the ocean. This buffering is achieved through a temporal smoothing of the freshwater flux (with a 10-year e-folding time) to avoid any spurious low-frequency variability in the freshwater input to the ocean. The flux is then integrated in three latitudinal bands (90–40°N, 40°N to 40°S, and 40–90°S) and passed to the ocean. In the north and in the tropical/subtropical bands, the flux is equally distributed over the ocean on the same latitudinal bands. In the south band, it is split in two contributions of 50% each corresponding to ice shelf melting and iceberg melting. The ice shelf melting is geographically and vertically distributed along Antarctica so as to mimic the observed distribution from Depoorter et al. (2013) as described in Mathiot et al. (2017). The iceberg melting is spread offshore following the observed geographical distribution of icebergs of Merino et al. (2016) and distributed vertically over the top 150 m, similarly to river runoffs.

Insufficient information is available from the atmospheric and land surface models on the temperatures of freshwater inputs to the ocean so a number of simplifying assumptions are made: The temperatures of rain and snow reaching the ocean are assumed to be that of the SST or the ice surface temperature in the ice-covered areas; the temperature of the riverflow is assumed to be that of the SST at the river mouth (except if the latter is negative, in which case riverflow is assumed to be at 0°C). The freshwater flux from iceberg melting is treated as runoff; hence, the latent heat required to melt the ice is ignored, and its temperature is set to the SST. In contrast, the freshwater flux from ice shelf melting is treated as ice at 0°C, and the latent heat required to melt it is accounted for.

The lack of representation of the energy content of precipitation, riverflow, and icebergs results in energy not being conserved exactly in the model. It should be noted that these are not the only non-energy-conserving processes in the model. A number of subgrid-scale parameterizations in the ocean (e.g., eddy-induced velocity, convection, horizontal momentum mixing, and turbulent kinetic energy dissipation), in the atmosphere (e.g., convection scheme), and at the ocean-atmosphere interface (e.g., wind stress interpolation) are not conserving energy exactly. A small lack of energy conservation is not a major issue as small energy sources and sinks do not prevent the model from equilibrating, at least if the nonconserving terms are stationary. Achieving a more exact energy conservation is an objective for the next version of our climate model.

Optional Model Components

Other model components can be activated in IPSL-CM6A-LR but are neither further described in this article nor used in the model experiments presented below. These include atmospheric chemistry/aerosol microphysics models such as the INteractions with Chemistry and Aerosols (INCA Hauglustaine et al., 2014), the REactive Processes Ruling the Ozone BUdget in the Stratosphere (REPROBUS Marchand et al., 2012), and the Sectional Stratospheric Sulfur Aerosol (S3A Kleinschmitt et al., 2017) models. Activation of one of these model components (instead of specifying atmospheric chemical composition and aerosol climatologies) requires a small retuning of either the LMDZ6A model or, in the case of S3A, the background stratospheric sulfur budget in order to ensure a similar baseline climate than in IPSL-CM6A-LR. The coupling of these chemistry and aerosol models with the other model components will be described in forthcoming publications.

Testing, Tuning, and Evaluation Procedure

The model was largely developed and tested under present-day climate using a setup we refer to as a pdControl setup, which corresponds to present-day climate forcings with an artificial sink of shortwave radiative energy reaching the ocean surface in order to compensate for the ongoing oceanic heat uptake of the current unequilibrated climate. A number of model features were tuned toward observations (see Hourdin et al., 2017, for the rationale).

The parameters considered during the tuning of the atmospheric model are given in Table 3 of Hourdin, Rio, Grandpeix, et al. (2020a). They concern in particular the control of the deep convection scheme, the control of the conversion of clouds condensed water to rainfall, and the control of the vertical dependency of the width of the subrid-scale water distribution for non convective clouds. A parameter was introduced as well in the “thermal plume” model to control the representation of the transition from cumulus to stratocumulus clouds (Hourdin et al., 2020a, 2020b). The threshold value for the conversion from liquid cloud water to rainfall as well as a parameter that controls the indirect effect of clouds were used for the final tuning of the global radiative balance because they affect specifically the optical thickness of liquid (low) clouds, thus modifying the total shortwave radiation much more than the longwave.

At some point during the development process, the main development stream switched from a pdControlto a piControl setup, which corresponds to preindustrial climate forcings. The spin up lasts several hundreds years but with some evolution of the model physics as the tuning was being finalized. The final tuning process involved changes in parameters associated with the sea ice physical properties (albedo and conductivity), subgrid-scale orography parametrization, and penetration of energy in the upper ocean with and without sea ice cover. Various options were envisaged as well concerning the control of atmospheric deep convection and its competition with shallow convection. One important choice of the final configuration was to consider boundary layer convective transport by the “thermal plume model” outside cold pools only. With this choice, thermal plumes see a more unstable environment (since the thermal plume is more stable than the mean column). Thermal plumes are therefore more active, which in turns favors shallow convection compared to deep. The atmospheric deep convection activity over the ocean was modified as well by using a different (larger) value of the horizontal density of cold pools: one cold pool per (33 km)² over ocean versus one per (350 km)² over land (note that the parameterization of this cold pool density is currently further tested in more recent versions). The surface drag over the ocean was also modified by introducing a gustiness term computed as a function of the vertical velocity associated with air lifting by the thermal plumes and by the gust fronts of cold pools.

The model code was then frozen (Version 6.1.0) and subsequently altered only for correcting diagnostics or allowing further options and configurations. Versions 6.1.0 to 6.1.11 (the current version) are therefore bit-reproducible for a given domain decomposition, compiling options, and supercomputer.

A multicentennial preindustrial control was then simulated: 100 years (1750–1850) as the piControl-spinup experiment and 2,000 years (1850–3849) as the piControl experiment. It should be noted that the piControl experiment suffers from a small cooling drift of ∼ 0.2 K in 2,000 years. A shorter piControl experiment of 250 years labeled r1i2p1f1 was run on the Joliot-Curie supercomputer to check the consistency. A large ensemble (32 members) of historical(1850–2014) simulations were performed following the CMIP6 protocol. Initial conditions for the historical simulations were sampled every 20 or 40 years of the piControl starting with Year 1870 of the piControl. The r1i1p1f1 simulation was selected qualitatively among the first ∼12 available members at the time of selection on the basis of a few key observables of the historical period such as the evolution of the global mean surface air temperature, summer sea ice extent in the Arctic ocean, and annual sea ice volume in the Arctic ocean. The rationale for highlighting a particular member is that we expect many users to only consider r1i1p1f1 rather than the whole ensemble. A more thorough ongoing analysis of our historical large ensemble shows that other members appear to be closer to the observed record in many respects. Most of the historical simulations were prolonged (outside the CMIP6 protocol) to 2059 using SSP245 atmospheric, land use, and solar forcings (except for the wood harvest and ozone forcings, not available at the time, which have been kept constant to their 2014 values). An ensemble of scenario simulations for 2015–2100 with a few extensions to the Year 2300 were also performed following ScenarioMIP guidelines (O'Neill et al., 2016).

Infrastructure Improvements

The IPSL-CM6A-LR model can be extracted, installed, and compiled on a specific machine using a suite of scripts called modipsl. Simulations are executed within the libIGCM running environment, which can be used to set up and run a simulation on a specific machine through a chain of computing and postprocessing jobs. Metadata from the simulations are sent to the Hermes supervising tool which can be used to monitor progress in the simulations, and key variables from different simulations can be intercompared using an intermonitoring tool on a dedicated web server (closed access).

The CMIP6 simulations were performed at the Très Grand Centre de Calcul (TGCC) on the Curie supercomputer with a switch during the CMIP6 production in October 2018 to the Joliot-Curie supercomputer. Both a piControl and a historical simulations initially performed on Curie were repeated on the Joliot-Curie supercomputer to ensure the climate statistics were comparable on both supercomputers. The throughput is about 13 and 16 simulated years per day on Curie and Joliot-Curie, respectively, on 960 processors with the full CMIP6 output.

Outputs from the IPSL-CM6A-LR model are managed by the XML Input/Output Server (XIOS Meurdesoif et al., 2016). For the CMIP6 production, model output variables in native format were kept to the minimum. Instead, and in sharp contrast to CMIP5, the CMIP6-compliant model output has been produced on the fly using XIOS methods from the code. XIOS is driven by XML files describing the whole netCDF file structure (dimensions, attributes, etc.). Such XML files were produced using the dr2XML (https://github.com/rigoudyg/dr2xml) python library developed by our CNRM-CERFACS collaborators, which translates the CMIP6 Data Request and Controlled Vocabulary (https://github.com/WCRP-CMIP/CMIP6_CVs) into XML files for XIOS. We use dr2XMLto generate XML files for each simulated year of a given CMIP6 experiment and member. The netCDF time series are created and filled by XIOS all along the simulation to avoid concatenation during postprocessing.

A quality assurance is applied at the end of each CMIP6 model simulation. The simulation is validated from a scientific point of view to make sure there is no critical issue or inconsistencies in the diagnostics (e.g., incorrect application of a forcing term, wrong sign, and recurring patterns). Each file then undergoes several checks against the CMIP6 controlled vocabulary to ensure its conformance with the CMIP6 Data Reference Syntax (https://pcmdi.llnl.gov/CMIP6) in terms of time axis and coverage, variable and global required metadata, filename syntax, etc.

The data are then published on the Earth System Grid Federation (ESGF), which guarantees a strong and effective data management. The esgprep toolbox (https://esgf.github.io/esgf-prepare) is a piece of software that eases data preparation according to CMIP6 conformance. Once a model simulation is validated and checked, the netCDF files are migrated in the proper CMIP6 directory structure with the esgprep commands in a shared space of the file system. The IPSL hosts an ESGF index with all data sets from the French climate simulations and a data node to disseminate data sets from the IPSL climate simulations. The IPSL CMIP6 data sets are published on the ESGF data node using the usual esgpublish command line provided by the node stack. During the publication process, the Persistent IDentifier (PID) included in each netCDF file, is permanently stored in a dedicated database at the German Climate Computing Centre (DKRZ) allowing further data citation.

Evaluation of Present-Day Climatology

In this section we evaluate the present-day climate of IPSL-CM6A-LR against our CMIP5 flagship configurations, IPSL-CM5A-LR and IPSL-CM5A-MR, by considering recent periods (i.e., the 1980–2005 period, when not mentioned otherwise) of our historical simulations. The references (i.e., observations and/or reanalyses) used to evaluate the models generally cover the same period, but may sometimes include some years after (like ERA-Interim) or before (like WOA13-v2). We argue that not considering the exact same periods for the simulations and the observations only has a minor impact on the results given that (i) the model internal variability is not synchronized with that of the observations and (ii) large volcanic eruptions are included in the periods considered. The list of evaluated model variables and data sets against which they are evaluated are presented in Table 1.

Table 1 List of Evaluated Model Variables and Data Sets Against Which They Are Evaluated

	CMIP6
Variable	variable		Period	Reference
full name	short name	Product	covered	publication
2-m air temperature	tas	ERA Interim	1980–2009	Dee et al. (2011)
3-D temperature	ta	—	—	—
3-D zonal wind	ua	—	—	—
3-D specific humidity	hus	—	—	—
3-D relative humidity	hur	—	—	—
Sea level pressure	slp or psl	—	—	—
10-m zonal and meridional wind component	uas, vas	—	—	—
Air temperature at 850 and 200 hPa	ta_850, ta_200	—	—	—
Zonal wind component at 850 and 200 hPa	ua_850, ua_200	—	—	—
Meridional wind component at 850 and 200 hPa	va_850, va_200	—	—	—
Geopotential height at 500 hPa	zg_500	—	—	—
Precipitation	pr	GPCP	1979–2009	Adler et al. (2018)
Precipitable water	prw	REMSS-PRW-v07r01	January 1988 to January 2019	Mears et al. (2018)
Longwave cloud radiative effect	rltcre	CERES-EBAF	2000–2012	Loeb et al. (2018)
Shortwave cloud radiative effect	rstcre	—	—	—
Upwelling shortwave at the top of atmosphere	rsut	—	—	—
Upwelling longwave at the top of atmosphere	rlut	—	—	—
Sea surface salinity	sos	WOA13-v2	1975–2004	Locarnini et al. (2013)
Atlantic meridional overturning stream function	msftyz	Smeed et al. (2017)	2004–2017	Smeed et al. (2017)
Northward oceanic heat transport	hfbasin	Ganachaud and Wunsch (2003)	1985–1996	Ganachaud and Wunsch (2003)

Thirty-two members have been performed for the historical period. Most of the diagnostics are not qualitatively sensitive to the choice of the member (when looking essentially at the mean state). We thus use only the first member (r1i1p1f1) in the diagnostics and illustrate the spread within the ensemble for some of the diagnostics. We present the most common variables used in climatology (like SST, surface air temperature, and precipitation) and concentrate on variables that are affected by the coupling between LMDZ, NEMO, and ORCHIDEE. Thus, we do not repeat the evaluation of variables that are close to those presented in LMDZ6A AMIP paper (Hourdin et al., 2020a). There have been numerous developments in the different components of the model; tracing back the evolution of the biases to particular developments requires a well-defined experimental framework (Bodas-Salcedo et al., 2019) and additional simulations which we have not performed. For this reason we focus this study on the evolution of the biases between the IPSL-CM5 and IPSL-CM6 models. The evaluation starts with surface temperatures (SST and surface air temperature) and follows with results for the atmosphere, the ocean and the sea ice.

SSTs and Surface Air Temperatures

We first evaluate SST simulated by the model keeping in mind that the average SST between 50°S and 50°N was tuned to fit the observations in the pdControl experiment. In this respect it should be noted that, toward the end of the development process, a slightly negative bias in the model SST was deliberately introduced, along with a tuning of some sea ice parameters, to partly compensate for a negative bias in summertime sea ice volume. The overall biases in the SST (Figure 1) have been significantly reduced between the IPSL-CM5A models and IPSL-CM6A-LR. Part of the improvement is due to the fact that IPSL-CM5A-LR was inadvertently tuned too cold. Nevertheless, the improvements are also clear when the mean bias is subtracted (figures not shown). The North Atlantic negative anomaly around 45°N associated with the position of the North Atlantic drift is slightly reduced with a value of −4.3°C in IPSL-CM6A-LR, compared to −6.8°C in IPSL-CM5A-LR (with the value taken as the minimum temperature from the 60–15°W, 40–55°N box). For comparison this index ranges from −6.8°C to +2.0°C (90% interval) with a median around −3.8°C in CMIP5 models, and ranges from −7.1°C to +1.1°C (90% interval) with a median around −3.7°C in CMIP6 models. We hypothesize that the increase in horizontal resolution (and the better representation of the ocean topography that comes with it) together with the improved atmospheric circulation in LMDZ (notably the influence of the orography) have contributed to improve the oceanic circulation and the resulting SST in the area. The East Boundary warm biases have also been reduced in both extent and amplitude in IPSL-CM6A-LR subsequent to the improvements in the boundary layer humidity and stratocumulus clouds in those areas (Hourdin et al., 2020a) and to a careful tuning of radiative fluxes. Yet this improvement is less clear in the tropical south Atlantic. In addition to global mean and latitudinal variation, the contrast between eastern tropical basins and the rest of the tropical oceans was used as a target, considering East Tropical Ocean Anomalies (ETOA) defined by Hourdin et al. (2015). A similar attention was given to the reduction of the classic latitudinal SST biases, which counteracts a tendency of the model to produce too cold midlatitude SSTs and a warm bias close to Antarctica (e.g., Wang et al., 2014). In the North Pacific, a warm bias (mostly during summer time) persists over the ocean in IPSL-CM6A-LR. This bias is much less visible in Figure 1 (top and middle panels) because as indicated above, the CMIP5 versions were globally too cold. Relative anomalies show that the North Pacific bias was already present, although slightly weaker. This bias, robust to many tests which were conducted during the tuning phase of the coupled model is also present in other CMIP5 and CMIP6 climate models (not shown). Its origin is still to be investigated. The cold bias over the equatorial Pacific, another classic deficiency of coupled models, is reinforced in IPSL-CM6A-LR.

Figure 1. Global distribution of the annually averaged sea surface temperature (SST) bias (in °C) for IPSL-CM5A-LR (upper panel), IPSL-CM5A-MR (middle panel), and IPSL-CM6A-LR (bottom panel). Biases are computed against data from the World Ocean Atlas (WOA13-v2 Locarnini et al., 2013). Global-mean biases, root-mean-square errors (RMSE), and correlation coefficients are provided for each model.

Nevertheless, the SST biases against observations are altogether significantly reduced, even when comparing to IPSL-CM5A-MR which uses the same atmospheric horizontal grid as IPSL-CM6A-LR, with a reduction of the root-mean-square error (RMSE) from 1.4 to 0.975 and an increase of the correlation coefficient from 0.986 to 0.988. It is difficult to assess the statistical significance of such subtle changes. However, it can be noted that the correlation for IPSL-CM5A-MR falls outside of the range from the IPSL-CM6A-LR ensemble members (0.9875 to 0.9885).

Consistent with the discussion above, there is a general reduction of the bias in surface air temperature (notably over the ocean) from IPSL-CM5A-LR to IPSL-CM6A-LR (see Figure 2). Globally the increase in resolution has surely played a role as can be seen by comparing IPSL-CM5A-LR and IPSL-CM6A-LR. However, the difference between IPSL-CM5A-MR and IPSL-CM6A-LR is largely attributable to improvements in the model physics and subsequent improvements in the radiative budget in LMDZ (Hourdin et al., 2020a) and a much more systematic and better tuning of the model key parameters in IPSL-CM6A-LR in order to adjust the radiative fluxes.

Figure 2. Distributions of biases in the surface (2-m) air temperature (in °C) for the annual mean (ANM; a, d, g), December-January-February (DJF; b, e, h), and June-July-August (JJA; c, f, i) for IPSL-CM5A-LR (a, b, c), IPSL-CM5A-MR (d, e, f), and IPSL-CM6A-LR (g, h, i). The bias maps are computed against ERA Interim reanalysis (Dee et al., 2011).

Compared with IPSL-CM5A-MR, the warm bias over the Amazon basin and tropical Africa is reduced in IPSL-CM6A-LR. Cheruy et al. (2020) attribute the improvement to the reduction of the overestimation of the SW downward radiation at the surface. The cold bias over Asia (especially in winter) is stronger in IPSL-CM6A-LR than in IPSL-CM5 models. The changes in the snow albedo in ORCHIDEE are likely to be the cause of this amplification of the bias in comparison to IPSL-CM5A-MR. The new snow scheme improves the realism of the physical properties of the snowpack (albedo, density) in IPSL-CM6A-LR relative to IPSL-CM5A-MR. However, due to strong surface-atmosphere couplings, the larger value of the snow albedo appears to cancel out a previous error compensation (Cheruy et al., 2020) and favors a too strong snow cover in these continental areas. Part of this deficiency in summertime may be explained by the fact that the parameterization of snow albedo does not account for shading effects in mountainous regions, a process which is thought to reduce the surface albedo on the scale of a model grid box. A strong negative bias is indeed observed over the Tibet including during summertime. A model development to account for the impact of orography on surface albedo is planned for a future model version.

In the northern high latitudes the biases have also largely changed, due primarily to the revision of the boundary layer scheme which allows more decoupling in stable situations. Modifications in the subgrid-scale orography parameters affecting the atmospheric circulation, the sea ice model, and the land surface scheme also contribute to the change. The warm bias over the northern part of Canada has been reduced in the annual mean. There is actually some compensation of a warm bias in summer and a cold bias in winter in IPSL-CM6A-LR that replaces a warm bias all year long in IPSL-CM5A-MR (Figure 2). The biases over the Arctic are linked to the position of the sea ice edge and depend to some extent on the member being considered. However, a large warm bias is consistently simulated in winter over the Arctic.

Over inland Antarctica, a cold bias can be seen in IPSL-CM6A-LR surface air temperature in all seasons. This cold bias is likely to correspond to a warm bias diagnosed in the reanalysis surface temperature from a comparison with weather station data (Fréville et al., 2014; Jones & Lister, 2015), but the magnitude of this cold bias (up to 8°C) exceeds the warm bias in ERA (up to 5°C). In LMDZ6, the boundary layer scheme was indeed improved to match the temperatures observed at Dome C (Vignon et al., 2018). LMDZ5 tended to prevent the decoupling of the surface from the atmosphere in very stable conditions (Cheruy et al., 2020). To reach a good agreement with the observations, the ice she et albedo was also changed in IPSL-CM6A-LR following Grenfell et al. (1994).

Consistent with the reduction in SW radiation bias, the strong warm summer bias in midlatitudes that was shared by many models participating in CMIP5 (Cheruy et al., 2014) is reduced in the CMIP6 version. However, it remains present in smaller areas, particularly on the Southern Great Plains. In these regions the bias results from complex interactions between the land surface and the atmosphere, especially through convection (Koster et al., 2004). It is also likely that the lack of parameterization of propagating mesoscale convective systems that are known to occur frequently in the region contributes to this bias (Moncrieff, 2019).

Atmospheric Variables

In this section we present the evaluation of a set of common atmospheric variables—namely, surface precipitation and wind, temperature, and atmospheric water on zonal mean diagnostics—and conclude with a set of evaluation metrics obtained with the PCMDI Metrics Package (PMP Gleckler et al., 2016).

Atmospheric Structure

The zonal mean temperature and zonal wind (Figure 3) show a decrease in the warm bias over the Antarctic and of the cold bias at 200 mb in the polar vortex at both poles. The cold bias at midlatitudes between 850 and 400 mb was present in IPSL-CM5A-LR, vanished in IPSL-CM5A-MR but reappears in IPSL-CM6A-LR. The most striking improvement is for the zonal atmospheric circulation (Figure 3, bottom row). The subtropical jets used to be too close to the equator in the two CMIP5 IPSL models. The difference in resolution between IPSL-CM5A-LR and IPSL-CM5A-MR led to only slight improvements. Despite the same horizontal resolution than IPSL-CM5A-MR, IPSL-CM6A-LR has a much better zonal circulation with jets moving poleward. This improvement is mainly due to the changes in the physics of the atmospheric model and the increase in vertical resolution from 39 to 79 layers.

Figure 3. Zonally and annually averaged air temperature (top row, in °C) and zonal wind component (bottom row, in m s−1) for IPSL-CM5A-LR (left column), IPSL-CM5A-MR (middle column), and IPSL-CM6A-LR (right column). The black contours show the ERA Interim reanalysis (Dee et al., 2011) climatology, and the model bias against ERA Interim is depicted by the color scale. The pressure on the vertical axis is expressed in Pa.

The atmosphere is more humid than in previous models (Figure 4): The specific humidity in IPSL-CM5A (both LR and MR) used to be too low (i.e., corresponding to a dry bias) in the lower troposphere in the tropics, and it is now slightly larger (i.e., corresponding to a wet bias) than in ERA Interim. In terms of relative humidity (RH), IPSL-CM6A-LR appears to be too saturated compared with ERA Interim between 30° and 60° in latitude (in both hemispheres). The wet RH bias in the free troposphere of the midlatitudes was already present to some extent in the previous versions. This bias is known to partly reduce with increasing horizontal resolution as illustrated by the comparison of the IPSL-CM5A-LR and IPSL-CM5A-MR versions, as well as the comparison between the LR and HighResMIP horizontal grid in stand-alone atmospheric simulations with the LMDZ6A version (Hourdin et al., 2020a). The main difference of the IPSL-CM6A-LR version compared to IPSL-CM5A-LR/MR is the much wetter lower troposphere, at around 800 hPa. This change is related to the parametrization of the boundary layer transport to the boundary layer top of the air evaporated at the surface which is much more efficient with the thermal plume model in the IPSL-CM6A-LR version than with the old eddy diffusion scheme. This contributed to dry the near-surface air over the ocean, in better agreement with observation, but also resulted in a moist bias in the lower troposphere when compared to ERA Interim. The subgrid-scale distribution of total (vapor and condensed) water within a grid box as a function of height may also play a role in this. This wet bias should be put in relation with the increase in equilibrium climate sensitivity (ECS) in IPSL-CM6A-LR relative to IPSL-CM5A-LR and the diagnosed increased contribution of the water vapor feedback to the ECS (see Section 6).

Figure 4. Zonally and annually averaged specific humidity (top row, in kg kg−1) and relative humidity (bottom row, in %) for the IPSL-CM5A-LR (left column), IPSL-CM5A-MR (middle column), and IPSL-CM6A-LR (right column). The black contours show the ERA Interim reanalysis (Dee et al., 2011) climatology, and the model bias against ERA Interim is depicted by the color scale. The pressure on the vertical axis is expressed in Pa.

Surface Precipitation

In terms of precipitation, biases are generally consistent between the three model versions (see Figure 5 for global maps), with the main changes concerning the tropics (see Figure 6).

Figure 5. Global distribution of the annual mean precipitation biases (in mm day−1) for the IPSL-CM5A-LR (upper panel), IPSL-CM5A-MR (middle panel), and IPSL-CM6A-LR (bottom panel) models. The bias maps are computed against the Global Precipitation Climatology Project (GPCP Adler et al., 2018).

Figure 6. Annual mean precipitation rate (in mm day−1) for (a) GPCP, (b) IPSL-CM5A-LR, (c) IPSL-CM5A-MR, and (d) IPSL-CM6A-LR in the tropical region. The climatology is computed over the 1980–2005 period for the models and 1980–2009 for GPCP (Adler et al., 2018).

The equatorial Pacific is dryer in IPSL-CM6A-LR, reinforcing a classic bias of coupled model, associated with the above mentioned negative SST bias. This dry bias, particularly strong over the Warm Pool, is probably one of the most negative aspect of this new model version. Preliminary analysis indicates that it may be associated with reduced surface evaporation as a consequence of the modification of boundary layer mixing by the thermal plume model in this region. In contrast, rainfall over the Maritime Continent is strongly overestimated. This strong overestimation, also present in stand-alone atmospheric simulations, seems to be related to parameters of the deep convection schemes, in particular those taking different values over ocean and over land, such as the vertical velocity at the basis of convective clouds and the density of cold pools. Preliminary analysis suggests that the improvement of the South Pacific Convergence Zone (SPCZ) seems to be related to the activation of the thermal plume model, and a better representation of the shallow versus the deep convective regimes. Meanwhile, the so-called double Intertropical Convergence Zone (ITCZ) issue, with overestimated rainfall south of the equator over the East Pacific, is less pronounced in the new version. The double ITCZ issue is sometimes associated to entrainment in convective clouds (Oueslati & Bellon, 2013, 2015). The rainfall is altogether reduced over the eastern part of tropical oceans due to the modification of the parameterization of stratocumulus clouds and a careful tuning of the parameters that control precipitation in these clouds.

Rainfall is generally increased over semiarid regions like north India, Sahel, Australia, or around the Mediterranean Sea, in better general agreement with observations. Another major improvement of the new version is the reduction of the strong dry biases over the Amazon basin (as stated above).

The global precipitation rate is overestimated in the last version of the model, more than in the previous versions. For CMIP5, the global rainfall was considered as a target of tuning, and a strong effort was done to reduce the mean rainfall, which otherwise was generally overestimated by the IPSL model, as is the case in most global climate models. This target was intentionally abandoned for the tuning of IPSL-CM6A-LR, which explains for a large part the overestimation by 0.3 mm day⁻¹ of the global precipitation rate (about 10% of the observed value). This positive bias in the global mean precipitation is common to many other models. It cannot be excluded that this overestimation is partly due to an underestimation of the observed precipitation rate attributable to an underestimation of light rain over the tropical oceans (Berg et al., 2010; Hourdin et al., 2020a; Stephens et al., 2012).

PMP Large-Scale Summary Statistics

This section provides a general synthetic view of the evolution of the climatology of the atmosphere of the IPSL models between CMIP5 and CMIP6. We have used the PCMDI Metrics Package (PMP) to calculate a set of large-scale performance metrics (Gleckler et al., 2008), also called summary statistics, to summarize the agreement between the climate simulated by the model over the recent period and a set of references (observations and reanalysis, as listed in Table 1). Figure 7 shows the results for the most common atmospheric variables: 2-m air temperature (tas), surface precipitation (pr), precipitable water (prw), pressure at sea level (psl), upwelling shortwave (rsut) and longwave (rlut) radiation at the top of the atmosphere, cloud radiative effect at top of atmosphere on longwave (rltcre) and shortwave (rstcre) radiation, temperature, zonal and meridional wind at 850 mb (ta850, ua850, and va850) and at 200 mb (ta200, ua200, and va200), and geopotential height at 500 mb (zg500).

Figure 7. Global metrics summarizing the performance of the IPSL-CM6A-LR members of the historical ensemble (blue lines) against IPSL-CM5A-LR (red line), IPSL-CM5A-MR (green line), and the CMIP5 multimodel ensemble (gray lines). The metrics are for the global monthly spatiotemporal root-mean-square error (upper panel), global annual mean bias (middle panel), and spatial correlation on the annual mean field (bottom panel) for 17 atmospheric variables (see full name of the variables in Table 1). The statistics are computed for the models on the 1980–2005 climatology against the reference data sets listed in Table 1. Note that the reference period for the observational data sets can be different from the period listed in Table 1, for example, 1989–2009 for ERA Interim, January 1979 to April 2018 for GPCP, and 2000–2018 for CERES-EBAF. Each model is represented by a line that connects the values of the metric obtained for the different model variables (vertical axes). For readability the columns are sorted so that the line connecting the IPSL-CM5A-MR results goes up from the left to the right. The metrics were computed with the PCMDI Metrics Package (PMP Gleckler et al., 2016).

We display the results of the metrics using parallel coordinates plots, which has the advantage to display raw results and avoid the necessary normalization of the portrait plot (Gleckler et al., 2008). For the sake of readability, the variables are sorted to display the results by increasing order of performance for the IPSL-CM5A-MR model. The individual members of IPSL-CM6A-LR (blue lines) are grouped together, with no single ensemble member coming out of the pack. For the large majority of the metrics, the results for IPSL-CM5A-LR (red line) and IPSL-CM5A-MR (green line) are out of the spread of the IPSL-CM6A-LR ensemble, demonstrating a difference in climatology that cannot be explained by internal variability. The RMSE computed over the globe over the 12 months of the climatological annual cycle (Figure 7, top panel) has decreased for all variables except for ta_850 and zg_500. For many variables (va_200, ua_200, psl, va_850, ua_850, rsut, rlut, rstcre and rltcre), the error has considerably decreased compared to IPSL-CM5A-MR and is within or at the bottom of the CMIP5 model range. The global bias has not necessarily decreased for all the variables. The colder atmosphere of IPSL-CM6A-LR compared with IPSL-CM5A-MR shown in Figure 3 explains the higher negative biases for ta_850 and ta_200. For tas, IPSL-CM6A-LR is a little colder than IPSL-CM5A-MR but still shows the benefits of a better tuning compared with IPSL-CM5A-LR (which was much colder). The global bias for the meridional wind at 200 mb (va_200) has also increased (it is more positive), when it is only slightly more negative at 850 mb (va_850), and much closer to 0 for the surface meridional wind (vas). For zg_500 the bias has increased (it is more negative) due to a general reduction of the altitude of the geopotential at this standard level, over the whole globe except the Antarctic (not shown). The global bias for the upwelling shortwave and longwave radiation at the top of the atmosphere has slightly increased in absolute value (actually close to IPSL-CM5A-LR). It has not really changed for prw, uas and va_850. The increase in the biases for ta_850 and zg_500 partly explains the relatively larger RMSE for these variables. The improvement is also striking when looking at the correlation coefficients (Figure 7, bottom panel) with all the variables experiencing higher correlations in terms of their annual mean patterns.

Oceanic Variables

This section evaluates the model in terms of sea surface salinity (SSS), the global vertical temperature profile, the structure of the Atlantic Meridional Overturning Circulation (AMOC) and mixed layer depth (MLD), the meridional heat transport and a set of mass transports through key transects (Figures 8–15).

Figure 8. Global distribution of the sea surface salinity (SSS, in %0) biases in annual mean for IPSL-CM5A-LR (upper panel), IPSL-CM5A-MR (middle panel), and IPSL-CM6A-LR (bottom panel). Biases are computed against data from the World Ocean Atlas (WOA13-v2 Locarnini et al., 2013).

Figure 9. Latitude-depth distribution of the global zonal mean ocean temperature bias (°C) for IPSL-CM5A-LR (upper panel), IPSL-CM5A-MR (middle panel), and IPSL-CM6A-LR (bottom panel). Biases are computed against data from the World Ocean Atlas (WOA13-v2 Locarnini et al., 2013) over the 1955–2015 period.

Figure 10. Annual monthly maximum of the mixed layer depth (MLD, in m) in the North Atlantic for the reconstruction of de Boyer Montégut et al. (2004; upper left panel), IPSL-CM5A-LR and IPSL-CM5A-MR (upper middle and right panels), and for three IPSL-CM6A-LR historical members: the first member (r1i1p1f1, lower left panel), a member visually identified with a shallow mixed layer (r8i1p1f1, lower middle panel), and a member visually identified with a deep mixed layer (r12i1p1f1, lower right panel).

Figure 11. Same as Figure 10 but for the Southern Ocean. For IPSL-CM6A-LR, the lower left panel is also r1i1p1f1, the lower middle panel is the member visually identified with a shallow mixed layer (r3i1p1f1), and the lower right panel is the member visually identified with a deep mixed layer (r17i1p1f1).

Figure 12. Vertical profile of the meridional overturning stream function (in Sv) at 26°N in the Atlantic Ocean for IPSL-CM6A-LR (red line), IPSL-CM5A-LR (dark blue line), IPSL-CM5A-MR (light blue line), and the RAPID-WATCH observations (black line Smeed et al., 2017).

Figure 13. Atlantic meridional overturning stream function (in Sv) for IPSL-CM5A-LR (top panel), IPSL-CM5A-MR (middle panel), and IPSL-CM6A-LR (bottom panel) as a function of depth and latitude, on average over the 1980–2005 period.

Figure 14. Barotropic stream function (in Sv) for IPSL-CM5A-LR (top panel), IPSL-CM5A-MR (middle panel), and IPSL-CM6A-LR (bottom panel).

Figure 15. (a) Global and (b) Atlantic Ocean meridional heat transport (in PW) for IPSL-CM5A-LR (dark blue line), IPSL-CM5A-MR (light blue line), and IPSL-CM6A-LR (red line) over the 1980–2005 period and corresponding direct observations (black stars with error bars) from Ganachaud and Wunsch (2003).

SSS

The climatology of SSS shows many evolutions since IPSL-CM5A-LR (see Figure 8). Overall, the SSS is globally reduced. This corresponds to a relative increase of the precipitation in subtropical basins. In the Atlantic Ocean, this translates into a reduction of the positive bias in the subtropics and an increase of the fresh bias in the subpolar latitudes. A SSS decrease is also visible in the South Pacific Ocean and Southern Ocean. The negative bias around Indonesia is corrected in IPSL-CM6A-LR in spite of an overestimation of precipitation locally. This may be due to enhanced exchanges between the Pacific and the Indian Oceans (see transport in Indonesian Throughflow in Table 2). The north and tropical parts of the Pacific Ocean are a little saltier in IPSL-CM6A-LR compared with IPSL-CM5A-MR, which is consistent with the reduction in precipitation in the area.

Table 2 Mass Transports (in Sv) Through a Selection of Key Transects of the Global Ocean as Defined in Griffies et al. (2016)

Transect	IPSL-CM5A-LR	IPSL-CM5A-MR	IPSL-CM6A-LR	Observations
Barents opening	−0.89	−0.059	4.06	2.0
Bering Strait	1.09	1.13	1.17	0.8
Denmark Strait	−5.47	−5.73	−5.26	−3.4±1.4
Drake Passage	101.46	109.2	150.87	136.7±6.9
Fram Strait	0.009	−0.86	−3.59	−2±2.7
Indonesian Throughflow	−10.72	−11.17	−13.60	−15
Mozambique Channel	−27.96	−27.21	−23.22	− 16.7±8.9

Note. In the three model configurations, the transports are computed as time averages over the period 1980–2005 of the historical simulations (r1i1p1f1 member). Mass transports are counted positively eastward and northward. Observations are from Griffies et al. (2016) and references therein.

Vertical Profile of Temperature

The vertical profile of temperature as a function of latitude has strongly evolved between IPSL-CM5A and IPSL-CM6A-LR (Figure 9). Overall, biases are larger in IPSL-CM6A-LR as compared to IPSL-CM5A, except to the north of 60°N, where warm anomalies are present in all versions. IPSL-CM6A-LR exhibit negative temperature anomalies in the Southern Ocean and globally below 1,500 m and positive ones above, except near the surface (see the discussion on SST in Section 3.1). Changes in the Southern Ocean are presumably associated with an increase in the ocean ventilation around Antarctica (Figure 11) and local negative surface air temperature anomalies in winter (Figure 2i). The cold ventilated water masses penetrate the deep ocean globally down to a depth of 2,000 m. Above, temperature anomalies are positive, reflecting the fact that the model is globally warmer in IPSL-CM6A-LR compared to IPSL-CM5A. Furthermore, the model presumably forms too much mode water, as found in many other climate models (Stouffer et al., 2017). A cold bias is also visible in the subtropical surface waters, reflecting the relatively cold SST (see Figure 1). Altogether, this can be interpreted as a stronger (weaker) thermocline (surface) stratification in midlatitudes in IPSL-CM6A-LR. To what extent this stronger stratification is an outcome of the tuning of the eORCA1 configuration used here, or a robust characteristic of the mean state in IPSL-CM6A-LR given the other components of the climate model, remains to be clarified.

Structure of the AMOC and MLD

We now turn to the oceanic general circulation and to the AMOC in particular. The above mentioned excessive thermocline stratification at midlatitudes translates into a pinching of the upper limb of the AMOC in the Atlantic Ocean (Figure 12). Indeed, at 26°N, the RAPID-WATCH observations suggest a maximum overturning around 1,000 m depth, while it is reached at 700 m depth in IPSL-CM6A-LR. In this respect, the vertical profile was more realistic in IPSL-CM5A configurations, but with a lower magnitude. Note that all versions of IPSL-CM exhibit an underestimation of the AMOC maximum at 26°N (by about 25% in IPSL-CM6A-LR), a bias that is common to many coarse resolution climate models in the absence of overflow parametrization (Danabasoglu et al., 2014). This may in part be explained by the difference in time period used in this comparison (2004–2017 for the observations as compared to 1980–2005 for the models). However, it is more likely to be due to biases in precipitation in the North Atlantic and/or the representation of overflows and western boundary currents, which remains a challenge in climate modeling.

The AMOC profile at 26°N also illustrates the excessive volume of cold deep water masses that is apparent in Figure 9: the stream function changes sign at a depth of around 2,800 m in IPSL-CM6A-LR versus 4,500 m in observations. Notwithstanding, the strength of the deep overturning cell is realistic, and the latitudinal extent of that cell compares well with previous versions of the model (Figure 13). In particular, the zero contour of the AMOC around 2,500 m depth, is very horizontal at all latitudes, a characteristic of all three model configurations. Above that contour, the positive AMOC cell, is maximum around 40°N in IPSL-CM6A-LR, as in previous versions of the model. This maximum reaches roughly 14 Sv in IPSL-CM6A-LR, which is notably larger than before. This may be related to the fact that dense water production in IPSL-CM6A-LR is different than in previous versions (Figure 10). In IPSL-CM5A, deep mixed layers are found south of Iceland and south of Greenland, which was unrealistic. This bias in IPSL-CM5A is associated to an over extended winter sea ice in the Labrador and Nordic Seas. In IPSL-CM6A-LR, deep mixed layers are confined to the Labrador Sea and the Nordic Seas, which is close to observed locations. Still, when looking at other members of the historical ensemble, it appears clearly that there is substantial variability in the North Atlantic deep convection in this model (Figure 10, the three panels of the bottom row).

Deep convection in the Nordic Seas may be directly related to the strengthening of the upper limb of the AMOC in IPSL-CM6A-LR north of 60°N (Figure 13, although the stream function at these latitudes were integrated along distorted model grid lines). The northward transport through the Barents opening is also more intense in IPSL-CM6A-LR compared to previous versions, and so is the return flow through the Fram Strait (Table 2). This can be interpreted as more intense exchanges between the North Atlantic and the Arctic, which is likely to affect sea ice there (see below).

In IPSL-CM6A-LR, deep convection in the Southern Hemisphere is also very intense, much more than in IPSL-CM5A models. There are important observational uncertainties related to the MLD estimations (Pellichero et al., 2017). However, this convection is possibly overestimated in IPSL-CM6A-LR (Figure  11). This provides cold water masses that invade the deep ocean and strengthens the meridional density gradients in the Southern Ocean, inducing a very strong Antarctic Circumpolar Current (Table 2, Drake Passage). This constitutes a major difference in barotropic stream function—and in the overall horizontal circulation—between IPSL-CM5A and IPSL-CM6A-LR (Figure 14).

One of the major influence of the ocean to the global climate is through the meridional heat transport. This quantity is closer to observations in IPSL-CM6A-LR compared to previous versions (Figure 15a), which is a substantial improvement. However, a strong convergence of heat at 40°S remains, a feature which was already present in IPSL-CM5A model versions but is likely to be unrealistic. Direct observations at that latitude are not available, but the common view is that the global meridional heat transport is southward in the whole Southern Hemisphere (Trenberth & Caron, 2001), which is not the case in our models. This seems to be related to the strong meridional gradient in density at that latitude, in particular in the Atlantic Ocean, a consequence of excessive mode water formation to the north, as described above. We also see an anomalous northward heat transport at 50°S in IPSL-CM6A-LR, presumably linked to the excessive dense water formation to the south.

In the Northern Hemisphere, the northward heat transport in IPSL-CM6A-LR is larger than in IPSL-CM5A versions. Notwithstanding, the simulated value remains slightly underestimated at the latitude where direct observations are available (24°N). Further north, it is very similar with observations, but this is due to an overestimated contribution from the Pacific Ocean (not shown). This might be partly responsible for the positive SST anomalies found in the north Pacific (Figure 1).

In the Atlantic Ocean, the meridional heat transport remains underestimated at all latitudes, particularly at tropical latitudes (Figure 15b). Still, this bias is much reduced in IPSL-CM6A-LR compared to IPSL-CM5A versions, which presumably contributes to reducing SAT biases over Europe and northern Africa (Figure 2).

Sea Ice

Arctic sea ice was one of the targets considered during the tuning process (Hourdin et al., 2020a). We targeted around 20,000 km³ of preindustrial annual mean Arctic sea ice volume and ultimately obtained slightly more, typically within a range 20,000–25,000 km³, much less than for IPSL-CM5A-LR, but a bit more than in IPSL-CM5A-MR. We also aimed for a seasonal cycle of ice coverage in our pdControl experiment that was broadly consistent with observations in both hemispheres. The Antarctic sea ice volume was not specifically considered during the tuning stage. Overall we obtain a reasonably realistic simulation of sea ice, significantly improved, as compared with IPSL-CM5A-LR.

In the Northern Hemisphere, there is less sea ice in IPSL-CM6A-LR than in both IPSL-CM5A models, which typically results in a better agreement with satellite data (Figure 16). The possible causes are a better tuning (Massonnet et al., 2018), higher model resolution, and more elaborated ice-ocean physics (Uotila et al., 2017; Vancoppenolle et al., 2009). Wintertime sea ice extent and area in IPSL-CM6A-LR slightly underestimate satellite retrievals, but is still within observational uncertainty. Regionally, there is a lack of winter sea ice in Okhotsk sea, associated with warm air temperatures, and less ice than observed in Barents Sea. Summertime area and extent are generally lower than observed, but are still within observational uncertainty. Excess summer ice decay occurs on the Siberian Shelf.

Figure 16. Time series of integrated sea ice diagnostics (area, extent—both in 106 km2—and volume—in 103 km3) over the Northern Hemisphere. (a) Sea ice area is the integral of ice fraction within a given region—here the Northern Hemisphere. (b) Ice extent is the total area enclosed within the 15% sea ice fraction contour. (c) Ice volume is the integral of ice fraction times thickness. IPSL-CM6A-LR simulations feature the historical r1i1p1f1 member in black, other historical ensemble members in gray (16–85% confidence interval), and selected scenario runs in color. The 16–85% confidence range is also shown for IPSL-CM5A-LR (blue) and IPSL-CM5A-MR (orange), for historical and RCP8.5 runs. Symbols depict passive microwave satellite-based retrievals from three different algorithms: Nasa Team (Cavalieri et al., 1996), Bootstrap (Comiso, 1996), and OSI-SAF (EUMETSAT Ocean and Sea Ice Satellite Application Facility, 1996). The upper and lower curves correspond to March and September, respectively.

Looking at both sea ice area and extent, the amplitude of the seasonal cycle appears to be on the high range. The annual mean volume and its seasonal cycle are within the rather wide observational range (Massonnet et al., 2018). There are noticeable simulated decadal fluctuations in sea ice volume.

In the Southern Hemisphere, IPSL-CM6A-LR overall improves over both previous CMIP5 models, in particular in summer (Figure 17). Wintertime sea ice extent is overestimated by 1–2 million square kilometers, sea ice area even more so. This points to the classic high concentration bias of current sea ice models. Summertime extent and area are within uncertainty range. As a result, the amplitude of the seasonal cycle of areal sea ice coverage appears to be somewhat overestimated. Sea ice volume varies between 5,000 to 25,000 km³ with the season in the preindustrial climate, which is much higher than in our CMIP5 models. It is mostly wintertime sea ice that decreases in the 21st century. Summertime sea ice also decreases, but less clearly.

Figure 17. Same as Figure 16 but for the Southern Hemisphere. The upper and lower curves correspond to the September and February months, respectively.

Model Evaluation From a CMIP5 Point of View

To complement the above evaluation of the model climatology, we now revisit the recommendations for CMIP6 made by Stouffer et al. (2017) based on the results of a survey made after the CMIP5 exercise. One of the main scientific challenges facing the climate modeling community (first reported in Meehl et al., 2014) is indeed to understand “[…] the origins and consequences of systematic model biases.” Stouffer et al. (2017) listed six main long-lasting (across the various CMIP exercises) model biases from the survey as major points for improvement: (1) the double ITCZ; (2) the Walker circulation, the dry Amazon basin bias, and tropical variability; (3) tropical and subtropical low clouds and the East Boundary warm bias; (4) a too deep tropical thermocline; (5) too warm and too dry continental surfaces during summertime; and (6) the position of the Southern Hemisphere subtropical jet. In this final subsection of the evaluation of the present-day climatology of the model we illustrate how some of these biases evolved between IPSL-CM5A and IPSL-CM6A-LR in order to focus the evaluation of the model on identified problems for the CMIP community.

We make use of the diagnostics described in the previous sections along with a set of large-scale evaluation metrics (Figures 18–21) presented in the context of other CMIP5 and CMIP6 climate models. Note that only CMIP6 models available on the ESGF at the time of writing this study have been considered. All model outputs were regridded to the same regular 3° × 2 ° resolution longitude-latitude grid before assessing global mean biases and RMSE against observations. The model ranking shows a consistent but varying improvement of IPSL-CM6A-LR over IPSL-CM5 for the metrics presented in the figures and discussed below. Some of these metrics have been considered during the tuning process of the model; hence, their improvement is expected. This is the case of the radiative metrics (OLR and OSR in Figure 18 and the SW and LW cloud radiative effects in Figure 19) even though the metrics used for the tuning are not exactly the same as those presented here. Indeed it was not the RMSE on the seasonal cycle (considered as metrics in the multimodel plots) which was used for tuning but rather its latitudinal dependency as well as the contrasts between East tropical oceans and the rest of the tropics (Hourdin et al., 2020a). Most of the tuning procedure aimed at reducing the main regional SST biases. Thus, the reduced RMSE on the SST is clearly an outcome of the tuning process.

Figure 18. Biases (crosses) and root-mean-square errors (RMSE, closed circles) against observations of SST averaged over the 65°S to 65°N latitudinal region (in °C), outgoing shortwave radiation (OSR, in W m−2), and outgoing longwave radiation (OLR, in W m−2) for CMIP5 models (in green) including IPSL-CM5A-LR and IPSL-CM5A-MR (in blue) and IPSL-CM5B-LR (in purple) and CMIP6 models (in black) including IPSL-CM6A-LR (in red). The averages of the scores for CMIP5 and CMIP6 models are shown as [less than]CMIP5> and [less than]CMIP6>, respectively. The models are ranked according to their RMSE. [less than]CMIP5> is shown as a range through random sampling of an equivalent number of CMIP6 models. The RMSE is computed on the mean seasonal cycle for the period 1979–2005 against SST from the input4MIPs data set and radiative fluxes from CERES-EBAF. In order to compare models with different native resolutions, fields are first interpolated on a regular 3° × 2° grid with a conservative regridding scheme before computing the global mean RMSE against observations.

Figure 19. Same as Figure 18 but for the shortwave cloud radiative effect (SW CRE, in W m−2), longwave cloud radiative effect (LW CRE, in W m−2), and precipitation (PR, in mm day−1). The RMSE and biases are computed against the radiative fluxes from CERES-EBAF and precipitation from GPCP.

Figure 20. Same as Figure 18 but for the double ITCZ index (DI Pacific, in mm day−1), Pacific precipitation (in mm day−1), and the Amazonian precipitation (in mm day−1). The DI Pacific index is the annually averaged precipitation over the southeastern Pacific Ocean in the [150–100°W, 20–0°S] region (this diagnostics is performed on the native grid). The Pacific precipitation is annually averaged over the [150°E to 100°W, 40°S to 40°N] region. The Amazonian precipitation is annually averaged over the [65–50°W, 15–0°S] region. The map shows the multimodel averaged CMIP6 precipitation bias, which guided the choice of the regions to compute the indices.

Figure 21. (left and middle panels) Mean June-July-August (JJA) near-surface temperature biases (in °C) for amip (circles) and historical (plus signs) CMIP simulations averaged over a North American [105–80°W, 35–50°N] region and a Eurasian [23–92°E, 45–58°N] region. The map shows the multimodel averaged CMIP6 JJA near-surface temperature bias, which guided the choice of the regions to compute the indices. (right panel) Same as Figure 18 but for the Southern Hemisphere jet latitude. The bias is computed from the annually and zonally averaged near-surface wind u10m‾ as a weighted average of the latitude ϕ over the [70°S, 10°S] latitudinal band: ∫70∘S10∘Smax(u10m‾−1,0)ϕdϕ/∫70∘S10∘Smax(u10m‾−1,0)dϕ. The RMSE on this graph is computed directly from u10m‾. The reference near-surface wind climatology is taken from the ERA-Interim reanalysis.

The rainfall and position of the jets were not directly considered as tuning targets because the results were seen as reasonable enough from the beginning, but if it would not have been the case, some additional work or tuning would probably have been done in this direction. Concerning the mean rainfall bias in Version IPSL-CM6A-LR, we already mentioned that it was abandoned as a target for tuning, explaining the increased bias compared to IPSL-CM5A and IPSL-CM5B. However, the bias of IPSL-CM6A-LR is only slightly larger than the averaged bias of CMIP6 models. The RMSE has slightly decreased as a result of the combination of regional decreases and increases in the errors, as discussed above.

Coming back to the six points listed by Stouffer et al. (2017), the following comments can be made:

• Progress made on the double ITCZ, although not a target for tuning, is illustrated on the maps of annual mean precipitation climatologies on Figure 6 and on the scores shown on Figure 20. It can be seen that the southern branch of the double ITCZ in the eastern part of the tropical Pacific basin (as well as in the tropical Atlantic) has weakened in IPSL-CM6A-LR compared to the IPSL-CM5A models. The scores on Figure 20 (left panel) show the improvement of the Double ITCZ Pacific Index, with IPSL-CM6A-LR (in red) getting closer to the observed value. It must be noticed however that this improvement is accompanied by a reinforcement of another classic bias over the Pacific Ocean, consisting in a cold and dry tongue over the equator that extends too far west toward the Maritime Continent.
• Concerning the Walker circulation, the dry Amazon basin bias, and tropical variability, we provide evidence for a reduction of the dry Amazon basin bias in Section 3.2.2 and Figure 20. We speculate that the improvement comes from a mix of better parameterizations of relevant local processes and a more realistic representation of the regional patterns of the radiative budget as teleconnections are known to influence precipitation in tropical South America (Yin et al., 2013).
• The East Boundary warm bias together with subtropical low cloud biases have also been reduced and are clearly one of the major improvements in IPSL-CM6A-LR (as noted in Section 3.1). This reduction is linked for a large part to the improvement of the representation of cumulus and stratocumulus clouds in LMDZ (Hourdin et al., 2019b) and to a careful tuning of radiative and latent heat fluxes at the surface over tropical oceans.
• Concerning the summertime warm bias over continents, IPSL-CM6A-LR shows almost no improvement against IPSL-CM5A-LR but a pronounced improvement against IPSL-CM5B-LR (see Figures 2c, 2f, and 2i). More importantly, the warm bias is much reduced in IPSL-CM6A-LR amip simulations for the two regions shown on Figure 21. It is thus clear that the lack of improvement between IPSL-CM5A-LR and IPSL-CM6A-LR is due to the general cold bias introduced by the tuning in the former version. The warm bias is reduced because of an improved—but still not perfect—shortwave radiative flux at the surface (Cheruy et al., 2014, 2020). Over the Southern Great Plains the complexity of land-atmosphere interactions together with the difficulty to represent the convective activity (Van Weverberg et al., 2018) in relation with the absence of representation of the propagating convection in the present model (Klein et al., 2006) can explain the remaining bias.
• The improvement to the position of the Southern Hemisphere subtropical jet is illustrated on Figure 21, with IPSL-CM6A-LR performing better than previous IPSL-CM5 model versions. The jets, which are located too close to the equator in most CMIP models, are known to generally move poleward when the horizontal resolution is increased (Hourdin et al., 2013a). This was clearly the case between the IPSL-CM5A-LR and IPSL-CM5A-MR model versions. However, IPSL-CM6A-LR—which has the same horizontal grid as IPSL-CM5A-MR—shows a much better location of those jets as seen in Figure 3 and, for the Southern Hemisphere, in Figure 21. We do not have a definite explanation so far, but it may be related to the much better tuning of the latitudinal dependency of radiative fluxes in IPSL-CM6A-LR, which in turn controls the thermal structure, itself tightly related to the zonal wind through the thermal wind balance.

Modes of Variability

We now turn to the main modes of variability of the model. We first present the El Niño–Southern Oscillation (ENSO) as it is the dominant coupled ocean-atmosphere mode of variability and the new behavior of the ocean multidecadal variability in IPSL-CM6A-LR and wintertime midlatitude variability and atmospheric blocking. We do not present the atmospheric modes of variability defined with the leading empirical orthogonal functions of dynamical variables like the North Atlantic Oscillation or the Pacific North America pattern because they are presented in a separate study together with the role of orography parameterizations on those modes.

ENSO

The ENSO is the leading mode of interannual climate variability, emerging from air-sea interactions in the tropical Pacific, but with climate impacts worldwide due to atmospheric teleconnections (e.g., Timmermann et al., 2018). In particular, decadal modulation of ENSO results in decadal fluctuations of the global mean surface temperature (GMST), a natural phenomenon that modulates anthropogenic climate change (e.g., Kosaka & Xie, 2013). It is thus a very important phenomenon to represent in a global climate model. Previous studies have underlined characteristic biases in the representation of the tropical Pacific climate, that translate into a misrepresentation of some key ENSO processes (e.g., Bayr et al., 2018). Such typical biases include a too strong equatorial upwelling (“the cold tongue bias”), excessively dry western equatorial Pacific, and the tendency to form a “double ITCZ.” The cold tongue and dry western Pacific biases have increased between IPSL-CM5A-LR and IPSL-CM6A-LR (Figure 1, top and bottom panels; Figure 5, top and bottom panels; and Figure 22) but the “double ITCZ” bias has been reduced, with a South Pacific Convergence Zone that extends less into the eastern Pacific (Figures 6b and 6d).

Figure 22. Longitude-time section of the Pacific Ocean 5°S to 5°N average mean seasonal cycle of SST (colors, °C), zonal wind stress ×1,000 (green contours, N m−2), and rainfall (black contours, mm day−1) for (a) observations, (b) IPSL-CM5A-LR, (c) IPSL-CM5A-MR, and (d) IPSL-CM6A-LR. Longitude-time section of the 5°S to 5°N typical anomalies during an ENSO event: SST (colors, °C °C−1), zonal wind stress × 1,000 (green contours, N m−2 °C−1), and rainfall (black contours, mm day−1 °C−1) for (e) observations, (f) IPSL-CM5A-LR, (g) IPSL-CM5A-MR, and (h) IPSL-CM6A-LR. The typical ENSO anomalies are obtained as a lead/lag regression to the November–January averaged Niño3.4 SST anomalies. The average of the 6 available IPSL-CM5A-LR, 3 available IPSL-CM5A-MR, and 32 available IPSL-CM6A-LR historical members is used. Observations (1980–2018) are from GPCPv2.3 (Adler et al., 2018) for rainfall and TropFlux (Praveen Kumar et al., 2012, 2013) for SST and zonal wind stress.

IPSL-CM5A-LR and IPSL-CM6A-LR reveal relatively similar ENSO pattern evolutions (Figures 22f and 22h), with events that tend to start too early in spring and display westward phase propagation unlike in observations, and end too late the following year (compare to Figure 22e). The cold and dry equatorial biases in the mean climate result in SST, wind and rainfall anomalies that are shifted west relative to those in observations, and too weak precipitation anomalies (e.g., Bayr et al., 2018). The amplitude of ENSO has increased by ∼ 40% between IPSL-CM5A-LR and IPSL-CM6A-LR (Figure 23), now being slightly above the observed value. One of the major issues of ENSO in IPSL-CM5A-LR was its seasonality (Bellenger et al., 2014). ENSO events indeed peak in boreal winter in observations, but IPSL-CM5A-LR tended to produce a maximum of equatorial Pacific SST variability in boreal spring (Figure 23), due to its tendency to produce events peaking in spring in addition to those captured on Figures 22e and 22f. This out of phase behavior has disappeared in IPSL-CM6A-LR, but the gap in amplitude between spring and winter ENSO signals remains too weak (Figure 23) and below the CMIP5 median.

Figure 23. Seasonally stratified standard deviation of the average Niño3.4 SST anomalies for observations (black), the 6 IPSL-CM5A-LR (dark blue), 3 IPSL-CM5A-MR (light blue), and 32 IPSL-CM6A-LR (red) members. The green shading indicates the range of values from the CMIP5 multimodel database (green curve = median, dark green = 25th to 75th percentiles, light green = 5th to 95th percentiles). Observations are from GPCP v2.3 (Adler et al., 2018) for rainfall and TropFlux (Praveen Kumar et al., 2012) for SST over the period 1980–2018.

Overall, some mean state biases thought to strongly influence ENSO representation have diminished (double ITCZ), while others have strengthened (cold tongue and dry equatorial biases). The major improvement in ENSO representation between IPSL-CM5A-LR and IPSL-CM6A-LR is a better representation of the ENSO seasonality, with other aspects of ENSO being quite similar in the two models.

Multidecadal Variability

The climate variability at decadal to multidecadal time scales has strongly evolved in the latest version of the model (Figure 24 showing 500 years of the piControl simulation of each of the model). The Atlantic Multidecadal Variability (AMV) index, defined as the time evolution of the SST anomaly averaged between 0 and 65°N in the North Atlantic, seems to be dominated by a longer time scale in IPSL-CM6A-LR compared to both IPSL-CM5A versions (Figure 24, top panel). IPSL-CM5A-LR is indeed characterized by a marked bidecadal variability (Escudier et al., 2013; Ortega et al., 2015), also present, yet with weaker intensity, in IPSL-CM5A-MR (Wen et al., 2016). In the new model, the typical AMV time scale is much longer: successive peaks in the AMV index are separated by about 200 years. This bicentennial variability is very robust to small modifications in the oceanic code (not shown) but weakens toward the end of our 1,200,year-long piControl simulation. It should be noted that a similar feature is found in at least another CMIP6 model (CNRM-CM6 Voldoire et al., 2019) that shares the same ocean model as IPSL-CM6A-LR. Exact origin of this behavior is still under investigation. The spatial pattern of the AMV (Figure 25) exhibits a strong subpolar center of action and a relatively weaker tropical one as compared to observations. Note however that the AMV pattern in HadISST (Figure 25) was computed from the 1920–2016 period and the global SST averaged between 60°S and 60°N was removed from all grid points before computing the North Atlantic average (0–60°N, 80–0°W) following (Trenberth & Shea, 2006). Hence, this represents variability over a shorter and different period than the 500 years of the preindustrial control and possibly still polluted by external forcings in spite of the detrending. The AMV pattern in IPSL-CM6A-LR is also marked by a relatively clear teleconnection in the Pacific, with a pattern resembling a negative phase of the Interdecadal Pacific Oscillation (IPO) associated with a positive AMV phase, as in observations. Both IPSL-CM5A models failed to reproduce this teleconnection.

Figure 24. Time evolution of the AMV index (top panel, in K), AMOC maximum (middle panel, in Sv), and Antarctic Circumpolar Current (ACC) measured as the mass transport through the Drake Passage (bottom panel, in Sv) in the IPSL-CM5A-LR, IPSL-CM5A-MR, and IPSL-CM6A-LR model versions. The AMV index is defined as the detrended 10-year low-pass filtered annual mean area-averaged SST anomalies over the North Atlantic basin (0–65°N, 80°W to 0°E). The AMOC maximum is taken from the meridional stream function between 10°N and 60°N and below 500 m. The mass transport at the Drake Passage is integrated from the surface to depth between the Cape Horn and the western Antarctic Peninsula.

Figure 25. Pattern of the Atlantic Multidecadal Variability (AMV) defined as the regression (in °C °C−1) of the global SST on the AMV index defined as the 10-year low-pass filtered annual mean area-averaged SST anomalies over the North Atlantic basin (0–60°N, 80°W to 0°E) for the HadISST observations (1920–2016 period) and the IPSL-CM5A-LR, IPSL-CM5A-MR, and IPSL-CM6A-LR models. The SST weighted average between 60°N and 60°S was subtracted from each grid point prior to any calculation, in order to account for the global warming trend following Trenberth and Shea (2006).

The difference in the main time scale of variability is also found in the evolution of the AMOC maximum (Figure 24, middle panel). In IPSL-CM6A-LR, the AMOC has a predominant variability at centennial time scales, with peak-to-peak amplitude of almost 4 Sv. The same index has weaker variability, and predominantly over a shorter time scale, in IPSL-CM5A versions. The intensity of the Antarctic Circumpolar Current (ACC) measured at the Drake Passage is also different between IPSL-CM6A-LR and IPSL-CM5A versions (Figure 24, bottom panel). In IPSL-CM6A-LR, there is a marked periodicity with an 80-year time scale, with peak-to-peak amplitude of up to 15 Sv. Such a periodicity is not visible in IPSL-CM5A models, although there seems to be a predominant variability at a similar time scale. The mechanisms leading to this variability are not yet fully understood. The AMOC centennial variability seems to be related to freshwater anomalies building up at very slow time scales in the Arctic Ocean and flushing into the North Atlantic Ocean. The links between AMV, AMOC, and ACC variability in IPSL-CM6A-LR remain to be investigated.

Wintertime Midlatitude Variability and Atmospheric Blocking

Figure 26a shows the frequency of wintertime blocked days in the r1i1p1f1 historical simulations of the IPSL models against observations. The envelope of blocking frequency from an ensemble of CMIP5 and CMIP6 models is also reported. Blocking is defined estimating the reversal of the daily geopotential height gradient at 500 hPa following D'Andrea et al. (1998). With respect to D'Andrea et al. (1998) here data is interpolated on a regular 2.5° × 2.5° grid, so that Δ=0°, ±2.5°, ±5°, and Φ_n=80°N, Φ₀=60°N, Φ_s=40°N. For a comprehensive review on blocking physics and climatology, the reader is referred to Woollings et al. (2018). It is particularly pertinent to analyze blocking frequency, because it has been a challenging phenomenon to reproduce for Numerical Weather Prediction (NWP) and global climate models alike for a long time. Davini and D'Andrea (2016) showed that there has been some improvement over generations of models, especially in the Pacific sector. In Europe, on the contrary, only a small number of models have blocking frequencies close to observed levels. This general tendency of climate models is by and large confirmed for the CMIP6 generation (see the light orange and light blue bands of Figure 26a). IPSL-CM6A-LR simulates more blocked days than the two IPSL-CM5 models over Europe (0–30°E) in better agreement with observations, although the frequency of blocked days is still underestimated. In this region, IPSL-CM6A-LR remains in line with the average behavior of the other CMIP6 models. There is a second maximum of blocking frequency at about 70°E, corresponding to Ural blocking, that is largely overestimated with respect to observations and other CMIP5 and CMIP6 models. The Pacific sector is also slightly overestimated.

Figure 26. Diagnostics of wintertime (DJFM) Northern Hemisphere midlatitude dynamics: (a) 1980–2005 frequency (in %) of blocked days as a function of longitude in observations (black line: average between JRA-55 Reanalysis, NCEP/NCAR Reanalysis, and ECMWF ERA-Interim Reanalysis), IPSL-CM5A-LR (dark blue line), IPSL-CM5A-MR (light blue line), and IPSL-CM6A-LR (red line). CMIP5 (shaded blue) and CMIP6 (shaded orange) multimodel ensemble spread is also shown as the ±1 standard deviation from the ensemble mean; (b) 1980–2005 mean zonal wind (in m s−1) at 500 hPa for IPSL-CM6A-LR; (c) stormtrack of IPSL-CM6A-LR, computed as the high frequency standard deviation (square high-pass filter at 6 days threshold) of the 500 hPa geopotential height (in m); (d) low-frequency standard deviation (6 days square low-pass filter) of the 500 hPa geopotential height (in m) for IPSL-CM6A-LR; (e) difference of mean zonal wind between IPSL-CM6A-LR and the ERA Interim reanalysis for the same period; (f) difference of IPSL-CM6A-LR stormtrack with ERA-Interim; (g) difference of IPSL-CM6A-LR low-frequency standard deviation with ERA-Interim. (h) difference of mean zonal wind between IPSL-CM6A-LR and the IPSL-CM5A-LR for the same period; (i) difference of IPSL-CM6A-LR stormtrack with IPSL-CM5A-LR; and (j) difference of IPSL-CM6A-LR low-frequency standard deviation with IPSL-CM5A-LR. Only the r1i1p1f1 historical simulations are considered.

In order to have a consistent understanding of the model behavior, Figures 26b–26g give an overview of the wintertime midlatitude variability of IPSL-CM6A-LR. Difference maps with IPSL-CM5A-LR are also shown in Figures 26h–26j. In the Atlantic sector, the midlatitude atmospheric jet is overestimated and too zonal, penetrating deeply into the European continent (Figures 26b and 26e) and carrying the Atlantic stormtrack along (Figures 26c and 26f). This brings about the underestimation of European blocking, and the overestimation of the Ural one. Over the Ural, excess low-frequency variability is consistently found (Figures 26d and 26g). The tendency toward excessively zonal midlatitude jets is linked to an underestimation of orographic drag (Pithan et al., 2016). In the Pacific sector the slight excess of blocking frequency is in agreement with a southward displacement of the jet (Figures 26b and 26e) and an excess cyclonic wave breaking (Rivière, 2009) at high latitudes, as visible in the variability maps (Figures 26d and 26g).

Improvements with respect to IPSL-CM5A-LR are clearly visible. The overestimation of the jet is much reduced in the new model (Figure 26h), which is consistent with the increase of blocking frequency in the Euro Atlantic sector. In IPSL-CM5 the jet is stronger and penetrates in the Eurasian continent slightly to the south with respect to IPSL-CM6A-LR. This causes larger low-frequency variability in IPSL-CM5 (Figure 26j) in a region spanning the eastern Mediterranean to the low latitudes of the the Siberian region. At the same time the southward displacement of the jet explains the absence of an overestimation Ural blocking frequency.

Simulations of the Historical Period Simulation of GMST

As a reminder, the members of our ensemble of historical simulations have the same natural and anthropogenic forcings and differ only in their initial conditions which were sampled every 20 to 40 years in the piControl simulation. All historical simulations have been prolonged to 2030 using SSP245 forcings. Because of the large uncertainties in the observations before the 1880s, the analysis here is limited to the 1880–2018 period.

Figure 27 shows the time evolution of GMST (here computed from the surface air temperature), both in absolute terms and as an anomaly relative to the 1880–2018 period. A large spread is present in both panels, with differences up to 0.75 K for a given year. The ensemble mean of the anomaly (Figure 27b) can be interpreted as the forced component of climate change (due to natural and anthropogenic forcings) with variations around it due to internal natural variability. We compare this anomaly to both the Cowtan and Way (2014) and (Rohde et al., 2013a, 2013b) observational data sets. The observed GMST time series is within the spread of the ensemble simulations but the model ensemble mean qualitatively departs from observed changes around 1935–1945 and since 2005 (Figure 27b). The departure from observations for the recent period is slightly enhanced if the anomaly is computed from the 1850–1899 reference period (not shown). This large range of possibilities in the GMST evolution of the historical members is induced by their different initial conditions. It is not restricted to interannual variability as the long-term warming trends also depends on the historical member. The observed GMST response to the Pinatubo volcanic eruption is well represented by the model ensemble mean but there are large differences in the GMST evolution in the period around the Pinatubo eruption (i.e., 1990–1994) depending on the phasing of natural modes of variability in the simulations.

Figure 27. Time evolution of (a) the annual global mean near-surface air temperature (GMST, in K) from the 31 historical members, prolonged using SSP245 forcings, of the IPSL-CM6A-LR model and (b) the anomalies of GMST relative to the 1880–2018 average for the Cowtan and Way (2014, in blue) and the Berkeley (Rohde et al., 2013a, 2013b, in green) data sets and for the model ensemble average (in black) of the individual historical members (in gray).

Figure 28 shows the observed and simulated recent warming trends over the 1978–2018 period. Some members compares better to observations, for example, with some degree of “warming hole” in the North Atlantic Ocean. There are however some discrepancies; in particular there is no member that reproduces the observed cooling trends in the Southeastern Pacific and the Southern Ocean to their full extent. The model generally reproduces the land/sea contrast in warming, with an average ratio of 1.61 (ranging from 1.52 to 1.79) between the global temperature over land and ocean over the 1978–2018 period compared to the observed ratio of 1.67 from the HadCRUT4 data set (Morice et al., 2013b) (the model data, SST over ocean and TAS over land, are regridded onto the observations temporally masked prior to the analysis). The Arctic amplification tends to be overestimated by the model, with an average trend of 0.88 K per decade over the 70–90°N region (ranging from 0.22 to 1.58 K per decade) relative to the global-mean trend of 0.26 K per decade (ranging from 0.16 to 0.36 K per decade), whereas the trend is about 0.79 K per decade over the 70–90°N region and about 0.19 K per decade for the global mean in the Cowtan and Way (2014) data set. Further work is going on to assess the diversity of historical members of the IPSL-CM6A-LR model and their relevance against observations.

Figure 28. Trends in near-surface air temperature (K yr−1) from the Cowtan and Way (2014) data set (upper left panel) and from the historical members of the IPSL-CM6A-LR model for the 1978–2018 period. Members are classified (from left to right, and top to bottom) by increasing spatial root-mean-square error (RMSE) relative to the observations. The member number and the RMSE values are indicated on the top left and bottom right corners of each panel, respectively. Dotted hatching indicates grid boxes where trends are significant (Mann-Kendall test, p[less than]0.1).

Carbon Fluxes

The historical simulations have prescribed CO₂ atmospheric mixing ratio as per observations (Meinshausen et al., 2017). Global fluxes to the ocean and land can be estimated from the spatially resolved flux calculations of the NEMO-PISCES and ORCHIDEE models in response of atmospheric CO₂ concentration and simulated climate (Figure 29 and Table 3). Compatible emissions are defined as the anthropogenic emissions that would be required to simulate the prescribed CO₂ concentration if the carbon cycle were to be fully interactive in the model. These compatible emissions can be diagnosed from the following equation: [Image Omitted. See PDF]where E_ff is the CO₂ emission flux from fossil fuel combustion and cement production, E_tot is the total anthropogenic CO₂ emission flux, E_lcc the CO₂ emission flux due to land cover changes (which is also estimated in the model), G_atm the growth rate of atmospheric CO₂ concentration, S_ocean the oceanic sink, and S_land the terrestrial sink (not accounting for changes in land cover).

Table 3 Decadal Mean Components of the Global CO₂ Budget for the 1990–1999 and 2009–2018 Periods for IPSL-CM6A-LR and From the Global Carbon Project (GCP Friedlingstein et al., 2019)

[Table omitted. See PDF]

Note.The GCP carbon budget shows an imbalance of 0.3 to 0.4 PgC yr⁻¹. The uncertainty ranges represent uncertainties for the GCP estimates and standard deviation across the 32 ensemble members for IPSL-CM6A-LR.

The ocean is a net sink of CO₂ and that sink increases from near 0 in 1850 to ∼2.9 PgC yr⁻¹ in 2018 with very little variability among the 32 historical members (standard deviation of ±0.07 PgC yr⁻¹). This simulated oceanic sink is consistent with the 2.5 ± 0.6 PgC yr⁻¹ estimate from the Global Carbon Project for the 2009–2018 decade (Friedlingstein et al., 2019). Similarly, the simulated oceanic sink over the 1990–1999 decade (2.1 ± 0.04 PgC yr⁻¹) is very similar to the 2.2 PgC yr⁻¹ flux diagnosed in IPSL-CM5A-LR.

The net terrestrial flux remains broadly negative (i.e., a source to the atmosphere) until approximately 1970 due to land cover change effects. The flux then increases and the land becomes a sink due primarily to the increasing CO₂ fertilization effect that dominates the land cover change effect. The net sink reaches 1.5 PgC yr⁻¹ in the last decade of the historical period. It should be noted that, over the 1990–1999 decade, the simulated sink is very close to that of the IPSL-CM5A-LR model (1.3±0.13 PgC yr⁻¹ compared to 1.28 ± 0.1 PgC yr⁻¹). The simulated net terrestrial sink is however less than the net flux of 2.1 ± 0.7 PgC yr⁻¹ estimated by (Friedlingstein et al., 2019) for the 2009–2018 decade. There is also a fairly large year-to-year variability due to climate variability at the regional scale (e.g., Schaefer et al., 2002), and a correspondingly large variability between the 32 ensemble members. The net terrestrial carbon fluxes, S_land−E_lcc, are consistent with the Global Carbon Project (Friedlingstein et al., 2019) estimates even though the E_lcc emissions are underestimated. This leads to simulated compatible emissions within the range of estimated fossil fuel emissions and cement production, E_ff, of the Global Carbon Project (Friedlingstein et al., 2019).

Figure 29. Prescribed atmospheric CO2 mixing ratio (ppmv, upper left), and inferred “compatible” emissions (upper right), oceanic (lower left), and terrestrial (lower right) net CO2 fluxes (in Pg C yr−1) in IPSL-CM6A-LR for the historical period (historical r1i1p1f1 in black and other ensemble members in gray) and the future scenario experiments (colored lines). A 5-year running average is applied for the fluxes. Estimates from the Global Carbon Project (Friedlingstein et al., 2019) and their uncertainties are indicated with black circles and gray error bars.

Transient Climate Response and ECS Estimates

Transient climate response (TCR) and ECS are two important quantities that characterize the model response to the CO₂ radiative forcing. There is increasing awareness however that these quantities are not intrinsic properties to the climate system (or to a given climate model) but may depend on the climate state (Mauritsen et al., 2019; Rugenstein et al., 2020). Furthermore estimates of these quantities depend on the details of how they are estimated in a particular model.

The ECS is traditionally defined as the equilibrium GMST change for a CO₂ doubling. We follow Gregory (2004) to estimate an effective ECS by assuming a linear—or quasi-linear—forcing-feedback relationship between the anomalies of the net downward radiative flux at the top of the atmosphere ΔR and the global mean surface air temperature ΔT to a giving forcing F, and extrapolating ΔT to its value for ΔR=F+λΔT=0, where λ is the feedback parameter.

To calculate these radiative and temperature anomalies, we subtract the preindustrial global mean value of the net downward radiative imbalance and near-surface air temperature (“rtmt” and “tas,” respectively, from the first 500 years of the preindustrial control run, 1850–2350) from the respective radiative and temperature values from the abrupt-2xCO2 and abrupt-4xCO2 r1i1p1f1 experiments. Results are sensitive to the length of the simulation because of the spatial and temporal dependence of the feedback parameters (Andrews et al., 2012, 2015; Knutti et al., 2017; Rugenstein et al., 2020). An ordinary least squares regression of the radiative imbalance on temperature anomalies results in equilibrium ΔT for CO₂ quadrupling of 9.05, 9.49, and 10.02 K depending on whether 150, 300, or 900 years of simulation are considered (Table  4). This translates into an ECS of 4.75 K (95% confidence interval: [4.32, 5.20]) when the fit is performed on 300 years and a factor of 2 is used to scale a quadrupling to a doubling CO₂ as per the usual assumption of a logarithmic dependence of the CO₂ radiative forcing on its atmospheric mixing ratio. This corresponds to a 17% increase from the value of 4.06 K in IPSL-CM5A-LR (95% confidence interval: [3.73, 4.41]), estimated by the same method (Table 4).

Table 4 Effective Radiative Forcing (ERF, in W m⁻²), Equilibrium Global Mean Surface Temperature Change (ΔT, in K), Different Estimates of the Equilibrium Climate Sensitivity (ECS, in K) as Derived From abrupt-2xCO2 and abrupt-4xCO2 Simulations Using Variants of the Gregory (2004) Method, and Transient Climate Response (TCR) for the IPSL-CM5A-LR and IPSL-CM6A-LR

Quantity/model	IPSL-CM5A-LR	IPSL-CM6A-LR
ERF 2×CO2 (W m−2)	—	3.50±0.27
ERF 4×CO2 (W m−2)	6.65±0.18	7.64±0.22
ΔT 4×CO2 (900 years, K)	—	10.02 (9.56, 10.62)
ΔT 4×CO2 (300 years, K)	8.12 (7.55, 8.74)	9.49 (8.65, 10.40)
ΔT 4×CO2 (150 years, K)	8.08 (7.36, 9.90)	9.05 (8.05, 10.20)
ECS from 4×CO2 (900 years, factor 2, K)	—	5.01 (4.76, 5.28)
ECS from 4×CO2 (300 years, factor 2, K)	4.06 (3.78, 4.37)	4.75 (4.33, 5.21)
ECS from 4×CO2 (150 years, factor 2, K)	4.04 (3.68, 4.45)	4.53 (4.02, 5.10)
ECS from 4×CO2 (300 years, scaled by ERF, K)	—	4.35
ECS from 4×CO2 (150 years, scaled by ERF, K)	—	4.15
ECS from 2×CO2 (300 years, K)	—	3.83 (3.03, 4.88)
TCR (K)	1.96 (2.09)	2.45

Note. ERF is calculated as in Lurton et al. (2020) by regressing the anomaly of the net radiative flux at the top of atmosphere against the anomaly in global mean surface temperature using the first 20 years of the experiment. The anomalies are computed after substracting the piControl values year by year. The confidence intervals correspond to ±2σ. For IPSL-CM5A-LR, we also provide for reference in parenthesis the TCR value published by Dufresne et al. (2013).

The factor of 2 scaling from 4 × CO₂ to 2 × CO₂ can be questioned because (i) the CO₂ forcing is supralogarithmic in the atmospheric CO₂ concentration as shown by Zhong and Haigh (2013) and Etminan et al. (2016) for detailed radiative transfer models and Lurton et al. (2020) for our climate model, and (ii) feedback parameters may depend on the forcing magnitude. In our model the effective radiative forcings (ERF) for CO₂ doubling and quadrupling are 3.46 and 7.53 W m⁻², respectively (Lurton et al., 2020). Renormalizing the surface temperature change extrapolated in the CO₂ quadrupling with the estimated ERF values leads to a reduced ECS for a CO₂ doubling of 4.35 K. This is to be compared to an ECS of 3.83 K if the regression is performed directly on the 300-year abrupt2xCO2 experiment. This last value may correspond better to the original definition of the ECS.

The transient climate response (TCR) is defined as the temperature change at the time of CO₂ doubling in an experiment where the CO₂ atmospheric concentration increases by 1% per year. More specifically it is computed as the 20-year global average in the 2-m surface temperature around the time of CO₂ doubling (years 61 to 80) in the 1pctCO2 experiment relative to the same quantity in the corresponding period of the piControl. TCR amounts to 2.45 K in IPSL-CM6A-LR against published values of 2.09, 2.05 and 1.52 K in IPSL-CM5A-LR, IPSL-CM5B-LR and IPSL-CM5B-MR, respectively (Dufresne et al., 2013). Thus, the larger ECS in IPSL-CM6A-LR also translates into a larger TCR in comparison to our previous generation of models.

Differences in ECS Between IPSL-CM5A-LR and IPSL-CM6A-LR

As discussed above, the effective ECS increases from 4.1 to 4.8 K between IPSL-CM5A-LR and IPSL-CM6A-LR. The relative contributions to ECS are calculated following Dufresne and Bony (2008) and Vial et al. (2013) and illustrated in the bar plots. This method decomposes the contributions to ECS into (i) rapid tropospheric and stratospheric adjustments to carbon dioxide and (ii) temperature-mediated feedbacks operating on longer time scales. More specifically the rapid tropospheric adjustment includes the climate response associated with all tropospheric adjustments (temperature, water vapor, and clouds), surface albedo change, and the small land surface warming due to the CO₂ forcing (Vial et al., 2013). The method also quantifies the relative contributions of the water vapor and temperature lapse rate, surface albedo, and cloud feedbacks. Individual feedbacks are calculated by the radiative kernel method (Bony et al., 2006; Soden et al., 2008; Shell et al., 2008). A radiative kernel acts as a partial derivative, representing the sensitivity of the radiative flux to changes in a climate variable, such as water vapor, temperature, and surface albedo. The radiative kernel is multiplied by the change in the climate variable of interest (i.e., water vapor) diagnosed from a model simulation and then normalized by the GMST change to produce the feedback value. We employ the same kernels as in Shell et al. (2008) for water vapor, temperature, and surface albedo. The cloud feedback is calculated as a corrected residual term, correcting for a cloud-masking term (Vial et al., 2013), which adds a consistent offset to net cloud feedback value estimated from the cloud radiative effect method (Andrews et al., 2012). A small residual term reflects nonlinearities in the relationship between radiative perturbation and the temperature response.

The main drivers of this larger ECS in IPSL-CM6A-LR are more positive rapid tropospheric adjustment to CO₂, and a stronger combined lapse rate and water vapor feedback (Figure 30a). We diagnose the strong tropospheric adjustment from aqua-4xCO2 and amip-4xCO2 simulations, as well as the abrupt-4xCO2simulations, and find that the stronger adjustments come from clear-sky regimes (not shown). The stronger water vapor feedback primarily results from strong moistening tendencies in weak ascent regimes around 500 hPa (Figure 30c). We diagnose this moistening tendency in weak ascent regimes by projecting the relative humidity anomalies, defined as the difference between relative humidity after 150 years of the abrupt-4xCO2 simulation and the piControl, into a circulation regime basis, wherein ω₅₀₀, the vertical pressure velocity at 500 hPa, acts as a proxy for the large-scale tropical circulation (Bony et al., 2004). This framework introduced by Bony et al. (2004) allows for attribution of changes in a climate variable to a given tropical circulation regime, ranging from strong ascent to strong subsidence regimes with increasing ω₅₀₀ values. Relative humidity anomalies reach up to 15% in these weak ascent regimes. However, it has also been shown that the IPSL-CM6A-LR model is too moist in the tropical atmosphere compared with ERA-Interim data (see Figure 4) which suggests this moistening might be exaggerated as well.

Figure 30. Diagrams supporting our analysis of the model's equilibrium climate sensitivity (ECS). (a) Bar plot showing the relative contributions (in K) to the ECS of the stratospheric adjustment, tropospheric rapid adjustments, combined lapse rate and water vapor (LR+WV), surface albedo, and cloud feedbacks for the IPSL-CM5A-LR and IPSL-CM6A-LR models. The residual term is due to nonlinearities in the feedback terms. (b) Anomalies in the tropical relative humidity (%) as a function of atmospheric pressure (hPa) and circulation regime as diagnosed by the vertical pressure velocity, ω500 in hPa day−1. (c) Distribution of the net cloud feedback (in W m−2 K−1) for IPSL-CM6A-LR. (d) Tropical net cloud feedback (in W m−2 K−1) as a function of the circulation regime as diagnosed by the vertical pressure velocity, ω500 in hPa day−1. The last two diagnostics are computed over the tropical ocean (30°N to 30°S).

The net cloud feedback, in contrast, is less positive in IPSL-CM6A-LR than in the previous model version. Compensating positive and negative feedbacks in the tropics give rise to a less positive tropical cloud feedback. Plotted in Figure 30b is the spatial distribution of the net global cloud feedback, calculated from the kernel method, in W m⁻² per K of GMST. A positive, warming feedback is in red, while a negative, stabilizing feedback is in blue. This feedback map demonstrates that the cloud response in IPSL-CM6A-LR is spatially heterogeneous, with large swathes of the tropical ocean covered in positive or negative cloud feedbacks. To interpret the spatial discontinuity between the regions of positive and negative cloud feedbacks, we project the net cloud feedback in the tropics onto the ω₅₀₀ basis, analogous to what was done for relative humidity anomalies. Based on the decomposition, the regions of positive cloud feedback can be linked to weak ascent regimes [−20, 0 hPa day⁻¹] and moderate to strong subsidence regimes [25, 100 hPa day⁻¹] (Figure 30d). By contrast, negative net cloud feedbacks arise in deep convective regimes and a portion of weak subsidence regimes. Moreover, we divide the net cloud feedback into SW and LW components to see whether the SW or LW component drives the net cloud feedback in particular regimes. In convecting regimes, the net cloud feedback more closely tracks the negative, LW cloud feedback, while in subsiding regimes, the net cloud feedback more closely follows the positive, SW cloud feedback (Figure 30d). The cloud feedback map shows that, geographically, positive values are found in regions of large-scale subsidence, which cover large parts of the tropical ocean and are associated with marine boundary layer cloud such as stratocumulus and shallow cumulus (Bony & Dufresne, 2005). By contrast, negative cloud feedback values occur in regions of deep convection, such as the Western Pacific Warm Pool. A negative feedback also occurs over the Southern Ocean, which could result from phase changes or thermodynamic changes with warming (Ceppi et al., 2016).

Future Scenarios Change in Surface Temperature

We now briefly present and discuss some results from the scenario simulations. The time evolution of the global mean surface air temperatures are shown on Figure 31. The temperature change in 2100 relative to 1850–1900 is larger than 2°C in all scenarios except the SSP119 where it overshoots 2°C before returning to below 2°C. It should be noted that the temperature change trajectory is very similar for all scenarios until circa 2040 when it starts to diverge according to the emission trajectory. This highlights the long time scales associated with the carbon cycle and the climate system (Collins et al., 2013). We also compare on Figure 31 the IPSL-CM5A-LR and IPSL-CM6A-LR models recognizing that the RCP and SSP scenarios are not fully equivalent as the repartition of the total net radiative forcing between the different terms has changed (Lurton et al., 2020). IPSL-CM6A-LR shows more warming than IPSL-CM5A-LR for the high-end scenarios (RCP245/SSP245, RCP6.0/SSP460, and RCP8.5/SSP585). This is expected from the larger TCR and ECS in IPSL-CM6A-LR. More surprising is the larger warming in IPSL-CM5A-LR for the historical period, which we attribute to a number of small differences in ERF. More specifically, the CO₂ ERF is smaller (1.59 vs. 1.83 W m⁻² in 2015), and on the contrary the ERF of the non-CO₂ greenhouse gases (CH₄, CFCs, N₂O, and O₃) is larger (1.58 vs. 1.03 W m⁻² in 2015) in IPSL-CM5A-LR compared to IPSL-CM6A-LR. The ERF for the anthropogenic aerosols is approximately the same for the two models (≈ −0.6 W m⁻²). Assuming that the climate feedback parameter and the ocean heat uptake efficiency are the same in the historical and 1pctCO2 experiments, we can indeed expect more warming in IPSL-CM5A-LR for the historical period compared to IPSL-CM6A-LR (1.61 vs. 1.43 K) despite a smaller TCR (2.09 vs. 2.45 K).

Figure 31. (upper panel) Change in global mean surface air temperature (GMST, in °C) relative to the 1850–1899 period in the historical r1i1p1f1 member (thick black line) and the other members (thin gray lines) and scenario experiments for the r1i1p1f1 member (thick colored lines) and other members (thin colored lines). Anomalies for 0°C, 1°C, and 2°C are indicated for reference. (lower panel) change in GMST relative to the 1850–1899 period for the IPSL-CM5A-LR and IPSL-CM6A-LR models for the historical period and the 21st century.

Distribution of Temperature and Precipitation Changes

Figure 32 shows the distributions of changes in surface air temperature normalized by the GMST change, for the SSP126 and SSP585 scenario experiments, both for the end of the 21st century (2070–2100 period), and at the end of the 23rd century (2270–2300, extended scenario runs). The normalized changes are defined relative to a 100-year preindustrial average. The patterns of change are quite similar for SSP126 and SSP585 at the end of the 21st century. In contrast, at the end of the 23rd century, patterns differ more between the two scenarios: SSP126 shows an Arctic warming pattern quite similar to that of 2100, whereas the relative warming for this region in the SSP585 scenario is less severe. However, SSP585 in 2300 shows an overall stronger warming in the Southern Hemisphere (if we average values on both hemispheres, we have a Northern Hemisphere to Southern Hemisphere ratio of 1.13:0.87 for SSP585 vs. 1.24:0.76 for SSP126), and its global patterns are more homogeneous than for the SSP126 scenario. The former shows more warming over the Southern Ocean while the latter exhibits a noticeable cold spot in the southern Pacific Ocean.

Figure 32. Geographical distributions of the normalized change in near-surface air temperature (in °C °C−1) for the CMIP6 SSP126 (upper panels) and SSP585 (lower panels) scenario experiments, at the end of the 21st century (2070–2100 period, left panels) and at the end of the 23rd century (2270–2300, right panels) as simulated by the IPSL-CM6A-LR model. The temperature change is defined relative to the preindustrial value (averaged over 100 years), and the normalization consists in dividing the local temperature change by the global mean surface air temperature change.

For precipitation, the patterns at the end of the 21st century are similar in both SSP126 and SSP585 scenarios, but they tend to differ more at the end of the 23rd century, with a somewhat smoother precipitation signature on some of the equatorial region for the SSP585 experiment (Figure 33).

Figure 33. Same as Figure 32 but for the relative change in surface precipitation normalized by the global mean surface temperature change (in % °C−1).

Changes in Sea Ice

We observe a rather large response of sea ice to the 21st century anthropogenic forcings, much larger than in IPSL-CM5A-LR. Summertime Arctic sea ice extent (Figure 16) responds more to changes in global mean temperature (SIMIP community, 2020). The simulated loss rate per °C of global mean temperature change in IPSL-CM6A-LR (−3.39 ± 0.87 × 10⁶ km² K⁻¹) has largely increased in comparison to IPSL-CM5A-LR (−1.48 ± 0.43 × 10⁶ km² K⁻¹) and IPSL-CM5A- MR (−1.67 ± 0.87 × 10⁶ km² K⁻¹). This is consistent with the near-zero summer Arctic sea ice extent for all scenarios in IPSL-CM6A-LR—a feature that is shared with the majority of CMIP6 models. It is also remarkable that winter sea ice almost disappears by 2100 in the fossil fuel intensive scenario (SSP585), which some of the other CMIP6 models also predict. Possible causes for this greater sensitivity, which should be further investigated, include the warm winter Arctic atmosphere, an ocean heat supply, changes in aerosol forcing, and ice drift. The ice volume loss starts in the early twentieth century and accelerates in its last three decades of the century. This is followed by a steady decrease over the 21st century. In the Southern Ocean (Figure 17), it is mostly winter sea ice that decreases in the 21st century. Summer sea ice also decreases, but less clearly so.

Changes in Carbon Fluxes

The land and oceanic net carbon fluxes for the scenarios are shown on Figure 29. The oceanic carbon uptake is projected to increase or decrease according to the scenario being considered, with a clear saturation occurring at large CO₂ concentrations (e.g., SSP585 and SSP460) and a decrease in the sink when the atmospheric CO₂ levels off or decreases. The net land carbon uptake peaks at about 3 PgC yr⁻¹ between 2020 and 2060 in all scenarios before returning to near 0, or even negative, values. The downregulation of the maximum photosynthetic capacity that was implemented to account for the impact of nutrient limitation on the CO₂ fertilization effect (see Section 4.2) may overestimate the limitation effect when atmospheric CO₂ concentration goes over 700 ppm (mainly after 2050) and thus may explain this extreme behavior at the end of the century. The formulation was chosen to broadly reproduce the change in gross primary production observed at Free Air Enrichment experiment when CO₂ is doubled (FACE Norby & Zak, 2011) but we overlooked the responses at very high CO₂ concentrations. A new parametrization is being designed and implemented as an option in the model.

Conclusions

We have described the main features of the IPSL-CM6A-LR climate model which has been developed at IPSL for CMIP6. We discuss the implementation of climate forcings in the model in Lurton et al. (2020) and will discuss the development philosophy and methodology in a future paper. In comparison to the previous generation of IPSL model several improvements have been introduced to the model: more physically based parameterizations (e.g., Hourdin, Rio, Grandpeix, et al., 2020a), more realistic implementation of some forcings (e.g., stratospheric aerosols), and more systematic tuning of adjustable parameters with a view to simulate key aspects of the model's climatology (SST, AMOC, and Arctic sea ice). The IPSL-CM6A-LR model performance is significantly improved over IPSL-CM5A-LR and IPSL-CM5A-MR and compares well to other published CMIP6 models for a number of metrics. However, some systematic regional biases and shortcomings persist (e.g., double ITCZ, frequency of midlatitude wintertime blockings, and ENSO dynamics).

The effective ECS (computed from a 300-year regression on abrupt-4xCO2 and divided by a factor of 2) increases from 4.1 to 4.8 K between IPSL-CM5A-LR and IPSL-CM6A-LR. The TCR correspondingly increases from 2.1 to 2.4 K. The increased ECS is due to increased contributions from tropospheric rapid adjustments and the combined lapse rate and water vapor feedback, which are only partly compensated by less positive cloud feedbacks.

A grand ensemble of 32 historical members has been performed with IPSL-CM6A-LR. The global mean surface air temperature increase simulated by the model is in the range 1.1 to 1.6 K in 2014 relative to 1850–1899 (across the ensemble members). While the ensemble mean warms more than the observations, some members are more consistent with observations. The IPSL-CM6A-LR shows a 1.6 to 6.8 K warming in 2100 across the scenarios relative to the same 1850–1899 period. The IPSL-CM6A-LR model exhibits a sea ice response to 21st century climate forcings on the high range in comparison to other CMIP5 and CMIP6 models.

A range of other papers in the Special Collection further evaluate particular aspects of the IPSL-CM6A-LR model. A comprehensive assessment of the model will require a lot more work in the coming years. We expect this to take place in the context of the CMIP6 multimodel ensemble on the basis of the vast amount of data we have published on the ESGF.

Conflict of Interest

The authors do not declare any competing interests.

Acknowledgments

The authors are grateful to the developers of the OASIS software and the Earth System model components not in the author list. This work was undertaken in the framework of the L-IPSL LABEX and the IPSL Climate Graduate School EUR. As such it benefited from the French state aid managed by the ANR under the “Investissements d'avenir” program with the reference ANR-11-IDEX-0004-17-EURE-0006. It also benefited from Belmont project GOTHAM, under Grant ANR-15-JCLI-0004-01, the ANR project ARISE under Grant ANR-18-CE01-0012, ANR project CONVERGENCE, under Grant ANR-13-MONU-0008 and MOPGA/Investissements d'Avenir project Archange, under Grant ANR-18-MPGA-0001. The CMIP6 project at IPSL used the HPC resources of TGCC under the allocations 2016-A0030107732, 2017-R0040110492, and 2018-R0040110492 (project gencmip6) provided by GENCI (Grand Équipement National de Calcul Intensif). This study benefited from the ESPRI (Ensemble de Services Pour la Recherche l'IPSL) computing and data center (https://mesocentre.ipsl.fr) which is supported by CNRS, Sorbonne Université, École Polytechnique, and CNES and through national and international grants. Support from the European Commission's Horizon 2020 Framework Programme is acknowledged, under Grant Agreement number 641816 for the “Coordinated Research in Earth Systems and Climate: Experiments, kNowledge, Dissemination and Outreach (CRESCENDO)” project (11/2015-10/2020) and under Grant Agreement number 820829 for the “Constraining uncertainty of multidecadal climate projections (CONSTRAIN)” project. Peter Gleckler and colleagues from the Program for Climate Model Diagnosis and Intercomparison (PCMDI) are acknowledged for their contribution to the performance metrics package. For the analyses we have used Python, CDO (https://code.mpimet.mpg.de/projects/cdo), NCL (http://ncl.ucar.edu/), R (https://www.R-project.org/) and took advantage of the CliMAF Python library (Climate Model Assessment Framework, https://github.com/rigoudyg/climaf).