Articles | Volume 15, issue 14
Model description paper
28 Jul 2022
Model description paper |  | 28 Jul 2022

The Earth system model CLIMBER-X v1.0 – Part 1: Climate model description and validation​​​​​​​​​​​​​​

Matteo Willeit, Andrey Ganopolski, Alexander Robinson, and Neil R. Edwards

The newly developed fast Earth system model CLIMBER-X is presented. The climate component of CLIMBER-X consists of a 2.5-D semi-empirical statistical–dynamical atmosphere model, a 3-D frictional–geostrophic ocean model, a dynamic–thermodynamic sea ice model and a land surface model. All the model components are discretized on a regular lat–long grid with a horizontal resolution of 5×5. The model has a throughput of  10 000 simulation years per day on a single node with 16 CPUs on a high-performance computer and is designed to simulate the evolution of the Earth system on temporal scales ranging from decades to >100 000 years. A comprehensive evaluation of the model performance for the present day and the historical period shows that CLIMBER-X is capable of realistically reproducing many observed climate characteristics, with results that generally lie within the range of state-of-the-art general circulation models. The analysis of model performance is complemented by a thorough assessment of climate feedbacks and model sensitivities to changes in external forcings and boundary conditions. Limitations and applicability of the model are critically discussed. CLIMBER-X also includes a detailed representation of the global carbon cycle and is coupled to an ice sheet model, which will be described in separate papers. CLIMBER-X is available as open-source code and is expected to be a useful tool for studying past climate changes and for the investigation of the long-term future evolution of the climate.

1 Introduction

Contemporary Earth system models (ESMs), based on relatively coarse-resolution (i.e. order of 200–300 km) atmosphere general circulation models (GCMs) and non-eddy resolving ocean models, currently reach simulation speeds of up to several hundred model years per day. This is sufficient to perform a single simulation or a small ensemble of simulations with the duration of several thousand years, the time typically needed to reach full equilibrium of the climate system and to perform near-future climate projections (IPCC, 2021​​​​​​​). However, there is also a need for much more computationally efficient climate models and ESMs suitable for performing simulations on orbital and even longer timescales as well for large multi-millennial ensembles of model simulations. This is why a certain ecological niche still remains for models of intermediate complexity, usually referred to as Earth system models of intermediate complexity (EMICs) (Claussen et al.2002). These computationally efficient models historically have rather coarse spatial resolution and, usually, employ simplified atmospheric components based on energy balance (i.e. UVic, Weaver et al.2001, BERN-3D, Ritz et al.2011), statistical–dynamical (CLIMBER-2, Petoukhov et al.2000, CLIMBER-3α, Montoya et al.2005), quasi-geostrophic (LOVECLIM, e.g. Goosse et al.2010) and 3-D (GENIE, Lenton et al.2007, and GENIE-PLASIM, Holden et al.2016) atmospheric models. Their oceanic components range from 2-D zonally averaged to full 3-D ocean general circulation models. Additionally, simplified and therefore fast GCMs have also been developed to fill the gap between EMICs and state-of-the-art GCMs (e.g. PlaSim, Fraedrich et al.2005, FAMOUS, Smith et al.2008, and ICCM, Farneti and Vallis2009). Since EMICs are positioned as Earth system models, they usually include a number of modules representing long-term processes, such as the global carbon cycle, ice sheets, or dynamical vegetation.

One such model is CLIMBER-2 (Ganopolski et al.2001; Petoukhov et al.2000). Its climate component was originally developed more than 20 years ago and has been widely used, primarily but not exclusively, for the study of past climates (e.g. Willeit et al.2019; Ganopolski and Brovkin2017; Ganopolski et al.2016). This model has a very coarse spatial resolution (minimal geographically explicit models), which had been dictated by the need to perform transient simulations on orbital timescales with the computers available at that time. Significant progress in computer performance made such coarse resolution unnecessary. At the same time, the availability of a vast number of both observational data and results of complex climate and Earth system models opens up the possibility of improving model performance and ensuring its consistency with more complex models.

Since on the one hand there is still demand for such models and, on the other hand, the availability of data and faster computers now enables the development of much better models but with a similarly low computational cost, we decided to develop a new model, CLIMBER-X. CLIMBER-X is based on a similar philosophy to CLIMBER-2 but with higher resolution, more internally consistent components and better treatment of individual processes.

A schematic illustration of all components of CLIMBER-X is presented in Fig. 1. Here we describe only the climate core of CLIMBER-X, which consists of the Semi-Empirical dynamical-Statistical Atmosphere Model (SESAM), the ocean model GOLDSTEIN, the SImple Sea Ice Model (SISIM) and the land surface–vegetation model PALADYN. CLIMBER-X also represents the global carbon cycle and therefore allows interactive simulation of the atmospheric CO2 and CH4 concentrations. Besides PALADYN, which simulates carbon cycle processes on land, the model includes HAMOCC (Maier-Reimer and Hasselmann1987; Ilyina et al.2013) as the ocean biogeochemistry and marine-sediment model. The ice sheet models SICOPOLIS (Greve1997) and Yelmo (Robinson et al.2020) are both included in CLIMBER-X as optional land ice components and are coupled with the climate component via an updated version of the physically based surface energy and mass balance interface SEMI (Calov et al.2005) and a basal ice shelf melt module. The viscoelastic mantle model VILMA (Klemann et al.2008; Martinec et al.2018) is used to simulate the bedrock response to changes in loading by solving the sea-level equation. The global carbon cycle model and ice sheet coupling will be described in detail in forthcoming papers.

2 Model description

The climate components of CLIMBER-X include SESAM, the frictional–geostrophic 3-D ocean model GOLDSTEIN, SISIM and the land model PALADYN (Fig. 1). They are coupled through fluxes of energy and water, the properties which are conserved in the entire system. Momentum is not conserved in the system, but surface wind stress is computed in the atmosphere and passed to the ocean component. The atmosphere, ocean, sea ice and land models are all discretized on the same regular longitude–latitude grid with a horizontal resolution of 5×5, which facilitates exchanges among them. The climate model components also share the same base time step of 1 d, which is also the coupling frequency between modules. Shorter time steps are used internally in the atmosphere and sea ice models for reasons of numerical stability. CLIMBER-X is designed to simulate the mean climatological state. It does not explicitly resolve the diurnal cycle and synoptic variability, although the effect of synoptic weather systems on the transport of energy and water is accounted for. The model does not represent interannual internal climate variability such as El Niño–Southern Oscillation (ENSO).

Figure 1Schematic illustration of the CLIMBER-X model, including exchanges and coupling between the different modules.


In CLIMBER-X the land–ocean distribution, topography and land cover can evolve with time. A high-resolution ( 10 arcmin) topography input file is used to automatically generate the land–sea mask, the topography, the orographic roughness and the ocean bathymetry. For the present day the RTopo-2 data of Schaffer et al. (2016) are used. Runoff routing is also derived automatically from the high-resolution topography by following the steepest slope and using the algorithm of Planchon and Darboux (2002) to fill depressions. No manual intervention is needed in the process, so that CLIMBER-X can in principle also be applied for simulations of times in the past when the distribution of continents was very different from the present day. CLIMBER-X also allows for dynamic transient changes in the land–sea mask, which is important for simulations on long timescales, when sea-level changes or ice sheet expansion and retreat can cause substantial changes to the land and ocean areas, with important effects on the fluxes of energy and water between the surface and the atmosphere. Furthermore, because the model resolution is relatively coarse and the removal or addition of a whole cell could cause a large perturbation in the model, we have implemented land–sea fractions in each grid cell to ensure more continuous transitions. This fractional approach, common to many models, also increases the effective horizontal model resolution.

Figure 2Meridional energy transport. (a) Total simulated meridional energy transport by ocean and atmosphere compared to observations (Trenberth and Caron2001; Fasullo and Trenberth2008). (b) Meridional energy transport by the atmosphere compared to estimates from Fasullo and Trenberth (2008) and modelled decomposition into dry static energy (DSE) and latent heat (LE) fluxes and contribution by eddies. (c) Global ocean meridional energy transport compared to observations (Fasullo and Trenberth2008).

CLIMBER-X is written entirely in Fortran 90+, making use of derived types to structure the single model components in a way which makes information transfer and exchange between different modules more transparent. The code is parallelized using the shared-memory OpenMP standard and integrates around 10 000 model years in 1 d on a single node with 2 × 8-core CPUs on a high-performance computer (Lenovo NeXtScale nx360M5, Xeon E5-2667v3 8C 3.2 GHz). The model is therefore  100–1000 times faster than state-of-the-art climate models that simulate weather when using comparable computational resources. Model input/output relies completely on the NetCDF data standard and uses the NCIO package (Robinson and Perrette2015) to read from and write to files. It therefore depends on the NetCDF library being available on the system where the model is run. Experimental settings and model parameters are defined in namelists, with one for each model component. A Python script is used to easily generate ensembles of simulations with different settings and parameters.

A conceptual overview of the individual model components is given in the following sections, while more details and equations are given in the Appendix.

2.1 Atmosphere model: SESAM

The atmospheric component of CLIMBER-X is designed as SESAM, based on a similar approach to the atmosphere of CLIMBER-2 (Petoukhov et al.2000). Originally, the atmospheric component of CLIMBER-2 was positioned as a model based on equations derived (with some assumptions) from first principles. However, it has been found that the governing equations and some key parameterizations are hard to derive from first principles, and a number of modifications have been necessary to match results with observational data or results of more physically based models. This approach is explicitly and extensively used during the development of SESAM. Compared to CLIMBER-2, the new atmospheric component of CLIMBER-X has a much higher spatial resolution (5×5), and most key parameterizations have been modified to improve performance against present-day climate and GCM sensitivity experiments. We made extensive use of available high-quality climatological data and results of complex climate model simulations which were not available when the CLIMBER-2 model was developed. This approach is reflected by the term “semi-empirical” in the model acronym.

Figure 3Annual mean (a) simulated meridional atmospheric streamfunction compared to (b) ERA-5 reanalysis (Hersbach et al.2020).

Figure 4Sea-level pressure in CLIMBER-X for (a) DJF and (d) JJA compared to ERA-Interim reanalysis (Dee et al.2011) (b, e). The zonal mean is additionally compared to CMIP5 models in panels (c) and (f).

The spatial resolution of SESAM, in terms of its representation of atmospheric fields, can be characterized as 2.5-D. This is related to the fact that, although prognostic variables (temperature, specific humidity, dust and eddy kinetic energy) as well as many diagnostic characteristics (cloudiness, condensation, shortwave radiation, etc.) are 2-D (latitude and longitude), the calculation of horizontal energy and water transport and the vertical fluxes of longwave radiation make use of information about the 3-D distribution of temperature, humidity, and wind velocity. All the vertical variations are purely diagnostic in the model. At the surface, each model grid cell is divided into different macro surface types, namely ocean, sea ice, land, lakes, and ice sheets. The macro surface types share the same atmospheric column but differ in surface elevation, albedo, and surface temperature and humidity. Sensible heat flux, evaporation, and shortwave and longwave radiation fluxes are computed separately for each macro surface type.

Figure 5Near-surface air temperature in CLIMBER-X for (a) DJF and (e) JJA and model bias relative to ERA-Interim reanalysis (Dee et al.2011) (b, f). The zonal mean is additionally compared to CMIP5 models in panels (c) and (g). The zonal mean absolute model bias in CLIMBER-X and a selection of CMIP5 models relative to ERA-Interim are shown in panels (d) and (h) for DJF and JJA, respectively.

The 3-D atmospheric structure is derived using assumptions about the universal vertical structure of temperature and relative humidity in the atmosphere. The temperature profile through the troposphere is a quadratic function of height, except in a 1500 m layer close to the surface where the lapse rate is a function of the near-surface atmospheric stability. Inversion strength is determined by the surface radiative balance and the difference between near-surface air temperature and skin temperature. Temperature is vertically uniform in the stratosphere. The tropopause height is derived assuming that the stratosphere is in radiative equilibrium. Relative humidity is taken to be constant within the planetary boundary layer, to decay exponentially in the low troposphere and to remain constant throughout the rest of the troposphere. Stratospheric relative humidity is prescribed at a uniform constant value. This approach differs from that used in CLIMBER-2, where specific humidity profiles were prescribed instead of relative humidity. This can lead to large biases in relative humidity in the mid–upper troposphere, where humidity is important for its effect on longwave radiation. Therefore, in SESAM the vertical profile of relative humidity is parameterized instead. See Appendix A1 for more details on the parameterization of the vertical structure.

Figure 6Zonal mean temperature simulated by CLIMBER-X for (a) DJF and (d) JJA compared to ERA-Interim reanalysis (Dee et al.2011) (b, e). The temperature bias relative to ERA-Interim is shown in panels (c) and (f). The dashed black line indicates the height of the tropopause.

Figure 7Near-surface air temperature seasonality in (a) CLIMBER-X compared to (b) ERA-Interim reanalysis (Dee et al.2011) computed as the difference between maximum and minimum monthly temperature.

Similarly to CLIMBER-2, in SESAM the horizontal wind velocity is divided into geostrophic and ageostrophic components. The geostrophic components of velocity at sea level are computed from sea-level pressure (SLP) and at any height through the troposphere by using the thermal wind approximation. The components of the surface wind are computed using the Taylor model (see e.g. Hansen et al.1983). The ageostrophic wind components in the planetary boundary layer are computed from the SLP gradient and the cross-isobar angle. SLP is the sum of zonally averaged and azonal components. The zonal component of SLP is derived based on the assumption that the strength and position of the major cells of the atmospheric meridional overturning circulation are controlled by average meridional temperature gradients and zonally averaged surface drag, similarly to Petoukhov et al. (2000). An additional dependence on zonal mean surface elevation has been added in SESAM to account for the effect of topography. The zonally averaged SLP is directly related to the zonal mean meridional (ageostrophic) wind in the planetary boundary layer and therefore the mean meridional circulation. The vertical structure of the zonally averaged meridional circulation is parameterized following the approach of Petoukhov et al. (2000). The azonal component of mean SLP resulting from stationary planetary waves is separated into a thermal component and an orographic component. The thermal component is assumed to be proportional to the azonal surface temperature anomaly reduced to sea level using the interrelation between long-term large-scale azonal temperature and pressure in quasi-stationary planetary-scale waves, similarly to Petoukhov et al. (2000). The effect of topographic stationary waves is particularly important for the Northern Hemisphere (NH) mid-latitude winter. The simple 1-D barotropic model for forced topographic Rossby waves of Charney and Eliassen (1949), applied to each latitudinal belt separately, is used to account for this effect in SESAM. Note that in CLIMBER-2 topographic waves were not represented because of the low spatial model resolution, whereas this limitation does not exist for CLIMBER-X. The main limitation in atmospheric dynamics in SESAM lies in the employed geostrophic approximation, which is not valid close to the Equator, resulting in deficiencies in the simulated dynamics in the tropics. See Appendix A2 for more details on the model dynamics.

SESAM includes a prognostic equation for vertically integrated energy and another for total-column water content, from which surface air temperature and humidity are derived. Sources and sinks of energy in the column include net longwave and shortwave radiation, surface sensible heat flux and latent heat release from condensation. Sources and sinks of water in the column are given by surface evaporation and precipitation. The 3-D wind field is used to advect the 3-D potential temperature field and 3-D specific humidity fields. Horizontal transport due to synoptic-scale processes is described as macroturbulent diffusion with an isotropic coefficient of horizontal diffusion expressed through the eddy kinetic energy (EKE). Precipitation is generated in the model whenever near-surface air relative humidity exceeds a critical threshold of 95 %. The water excess is then instantaneously added to precipitation. Over land, precipitation is additionally generated from column water content with a specified turnover time that is inversely proportional to atmospheric relative humidity. See Appendices A3 and A4 for more details on the energy and water balance equations.

Figure 8Zonal mean relative humidity simulated by CLIMBER-X for (a) DJF and (c) JJA compared to ERA-Interim reanalysis (Dee et al.2011) (b, d). The dashed black line indicates the height of the tropopause.

Figure 9Precipitation modelled by CLIMBER-X for (a) DJF and (d) JJA compared to ERA-Interim reanalysis (Dee et al.2011) (b, e). The zonal mean is additionally compared to CMIP5 models in panels (c) and (f).

Figure 10Simulated zonal mean cloud fraction and top-of-the-atmosphere net radiation fluxes for (top) DJF and (bottom) JJA compared to satellite observations (Loeb et al.2018), reanalysis (Dee et al.2011), and CMIP5 models.

The (vertically integrated) EKE is a prognostic variable in the model. EKE production is proportional to the baroclinicity of the atmosphere as defined by the Eady baroclinicity measure (e.g. Hoskins and Valdes1990). EKE dissipation depends on the surface aerodynamic drag coefficient and is proportional to EKE3/2. EKE is itself transported by macroturbulent diffusion and advected by the wind at 700 hPa. The synoptic component of surface wind speed is then taken to be proportional to EKE. See Appendix A5 for more details on the model representation of synoptic processes.

SESAM includes a single effective cloud layer. In this respect it is therefore even simpler than CLIMBER-2, which separately treats stratiform and convective clouds. We did not find any advantages in a separate treatment of different types of clouds in the framework of our modelling approach. Cloud fraction is a function of atmospheric relative humidity, effective vertical velocity and near-surface temperature inversion strength. Cloud top height is related to the height of the troposphere, modified by a factor depending on the vertical velocity at the cloud base. Cloud optical thickness is a function of surface air temperature, cloud fraction, column water content and sulfate aerosol load. Since CLIMBER-X does not include an atmospheric chemistry module, the spatial distribution of sulfate aerosol load must be prescribed as input to the model. See Appendix A6 for more details on cloud parameterizations.

The shortwave radiation scheme accounts for water vapour, clouds, sulfate aerosols, dust and ozone. The total shortwave radiation range is divided into two subintervals: ultraviolet + visible and near-infrared. The cloud albedo is a function of solar zenith angle and the optical thickness of the clouds (Feigelson et al.1975). The optical properties of dust are treated following the Yamamoto and Tanaka (1972) scheme. The integral transmission function of ozone is taken from Lacis and Hansen (1974). In the near-infrared band the absorption due to water vapour is described according to Feigelson et al. (1975). The effect of sulfate aerosols is included following the approach of Bauer et al. (2008). The scheme accounts for both the direct radiative forcing and the indirect effect through changes in cloud optical thickness. For simplicity, the radiative effect of black carbon and organic carbon aerosols is not represented in the model, as it is implicitly assumed that their combined net effect is negligible to a first approximation. See Appendix A7 for more details on shortwave radiation.

Figure 11Simulated ocean overturning circulation: (a) global and (b) Atlantic.


Figure 12Vertical profile of the simulated Atlantic meridional overturning streamfunction at 26 N (black) compared to observations from the RAPID array (Frajka-Williams et al.2021) (blue) and a selection of CMIP6 models (grey). The CLIMBER-X and CMIP6 streamfunction is computed from historical simulations as the average over the time period from 2000 to 2014, while the RAPID values represent an average from 2004 to 2020.

Figure 13Maximum monthly mixed-layer depth simulated by CLIMBER-X (a) compared to ECCOv4 reanalysis (ECCO Consortium et al.2021) (b).

To compute the longwave radiation fluxes, the atmosphere column is subdivided into 15 unevenly spaced levels. The longwave radiative scheme accounts explicitly for water vapour, carbon dioxide, and ozone through their integral transmission functions. The longwave radiative effect of the well-mixed greenhouse gases CH4, N2O and CFCs is accounted for through the use of an equivalent CO2 concentration in the radiative scheme, using the respective radiative forcing following Etminan et al. (2016). Similarly to sulfate aerosols, the 3-D field of ozone concentration must be prescribed as input to the model. It can be time-dependent. See Appendix A8 for more details on longwave radiation.

A simple representation of the global dust cycle is included in SESAM based on a single representative dust particle diameter (Bauer and Ganopolski2010). Dust emissions are computed in the land model and are an input to the atmosphere, where the dust concentration is assumed to decrease exponentially with height with a fixed height scale of 2 km. Dust is transported by the 3-D wind field and by macroturbulent diffusion. Dust sinks include both dry and wet deposition.

Figure 14Zonal mean ocean temperature simulated by CLIMBER-X (a, d, g) compared to WOA13 data (Levitus et al.2015) (b, e, h) for the (a, b, c) Atlantic, (d, e, f) Pacific, and (g, h, i) Indian and Southern oceans. Panels (c), (f), and (i) show the model bias.​​​​​​​

Figure 15Zonal mean ocean salinity simulated by CLIMBER-X (a, d, g) compared to WOA13 data (Levitus et al.2015) (b, e, h) for the (a, b, c) Atlantic, (d, e, f) Pacific, and (g, h, i) Indian and Southern oceans. The right panels show the model bias.

2.2 Ocean model: GOLDSTEIN

The ocean component of CLIMBER-X is based on GOLDSTEIN (Edwards et al.1998; Edwards and Shepherd2002; Edwards and Marsh2005). GOLDSTEIN is a 3-D frictional–geostrophic balance model and is also employed in several other existing EMICs, e.g. GENIE (Marsh et al.2011) and Bern3D (Müller et al.2006). In the original GOLDSTEIN implementation, the equations were coded in terms of non-dimensional quantities. For the inclusion in CLIMBER-X, the code has been dimensionalized for better readability and more efficient debugging.

The horizontal velocity in the ocean is diagnosed from a frictional–geostrophic balance and the vertical velocity is then derived from the continuity equation. The model assumes hydrostatic balance. The velocity relaxation to the velocities of the previous time step employed in GENIE and Bern3D has been removed in CLIMBER-X without drawbacks for model stability. Islands for barotropic flow are automatically derived from the topography, but the actual islands around which barotropic flow is permitted can be controlled by namelist settings. By default, non-zero barotropic flow is allowed through the Drake Passage and through the Indonesian archipelago. The Bering and Davis straits are closed for barotropic flow, but baroclinic tracer exchange is permitted. The friction for baroclinic velocities is set to be 4 d−1 and is globally uniform, while the friction for barotropic velocities is increased by a factor of 3 near continental boundaries and shallow topographic features. Close to the Equator the minimum absolute value of the Coriolis parameter is set to 5×10-6s−1. In contrast to other models using GOLDSTEIN (e.g. Holden et al.2016; Marsh et al.2011; Müller et al.2006) in CLIMBER-X, the wind stress as simulated by the atmospheric model is applied at the ocean surface, without any enhancement factor. Seawater density is computed using the UNESCO equation of state of Millero and Poisson (1981).

Figure 16Seasonal variation of total sea ice area for (a) the NH and (b) the SH as simulated by CLIMBER-X (black) compared to observations (Meier et al.2021) (blue) and CMIP6 models (grey).

Figure 17Sea ice concentration in the NH and SH in CLIMBER-X (top) compared to observations (Meier et al.2021) (bottom) for (left) March and (right) September.

The principal limitations to the ocean dynamics follow from reduced resolution and the neglect of non-linear momentum advection. The result is that there are no ocean eddies and hence no mechanism for internally generated variability and that boundary currents are unrealistically broad. One notable consequence of the strong momentum damping is a generally weak Antarctic Circumpolar Current and Gulf Stream. See Edwards and Marsh (2005) and Müller et al. (2006) for a more detailed discussion of the limitations of the frictional–geostrophic dynamics.

Figure 18Modelled permafrost extent and active layer thickness compared to the observed extent of continuous, discontinuous, and isolated permafrost (red lines, from dark red to light red) from Brown et al. (1998). The active layer thickness is calculated as the mean over the period 1981–2010 in grid cells that are permafrost during the whole time period.

The tracer transport is described by an advection–diffusion equation, written in flux form. To reduce numerical diffusion associated with the original advection scheme, we have implemented a flux-corrected transport scheme following Zalesak (1979). The model uses an explicit isopycnal and diapycnal diffusion scheme in its small-angle approximation (Redi1982) combined with a diffusive parameterization of the eddy-induced transport following the Gent–McWilliams skew flux (Gent and Mcwilliams1990; Griffies1998). The coefficients for both eddy-induced transport and isopycnal mixing are assumed to be equal by default (1500 m2s−1). Slopes of the isoneutral surfaces above 1×10-3 are tapered following Gerdes et al. (1991). The diapycnal diffusivity profile follows Bryan and Lewis (1979) with a minimum diffusivity of 0.1 cm2 s−1 at the surface and a maximum of 1.5 cm2 s−1 at the ocean floor. If the stratification of a water column is statically unstable, convective adjustment is applied using the scheme of Rahmstorf (1993). A simple mixed-layer scheme, based on Kraus and Turner (1967), is used in the model. Heat and freshwater fluxes at the ocean surface provide the top boundary conditions for the temperature and salinity prognostic equations. Since GOLDSTEIN is based on the rigid-lid approximation, the freshwater flux is implemented in terms of a virtual salinity flux, computed using local salinity. A global correction is employed to ensure that the global annual net freshwater flux into the ocean is zero, in the absence of changes in land ice volume. This is not sufficient to ensure conservation of salinity in the ocean because of the anti-correlation between surface freshwater flux and surface salinity (e.g. Yin et al.2010). Hence, the virtual salinity flux is additionally corrected to ensure that the annual mean global net surface flux is zero. The geothermal heat flux of Lucazeau (2019) is applied as a bottom ocean boundary condition.

The model is discretized on a regular lat–long grid with 5×5 resolution, matching the atmosphere model resolution, and 23 unequally spaced vertical layers, with a 10 m top layer and layer thickness increasing with depth and reaching 500 m at the ocean bottom. Ocean bathymetry is smoothed by convolving bathymetry with a directly neighbouring four-grid-point kernel to remove strong gradients that would lead to numerical instabilities. To improve model stability, an additional limitation on the topographic slope entering in the Coriolis term (planetary vorticity advection) of the barotropic velocity equation is applied.

Since the land–sea mask is fully flexible in CLIMBER-X, GOLDSTEIN has been adapted to deal with newly forming and disappearing oceanic grid cells. It is always ensured that all grid cells of the ocean domain are connected. If “ocean lakes” are formed at any point in time, the isolated ocean cells are removed from the ocean model domain and are treated as lake in the land model instead. A minimum ocean fraction of 0.1 in a 5×5 grid cell is required for a cell to be considered part of the ocean domain. Newly formed cells are initialized using information from neighbouring grid cells. Changes in sea level, and therefore ocean volume, are additionally accounted for by scaling the thicknesses of the ocean layers below a depth of 1000 m to match the actual ocean volume derived from the high-resolution topography and provided as input to the ocean model. Total tracer inventories in the ocean are conserved in this process. Islands for barotropic flow are automatically updated with the mask. Velocities are diagnostic quantities in the model, and therefore no intervention on velocity fields is required when updating the land–sea mask.

Information on sea surface height is needed in the sea ice momentum equation. The rigid-lid formulation employed in GOLDSTEIN does imply a surface pressure (e.g. Pinardi et al.1995), the so-called rigid-lid pressure, which is directly related to the sea surface height. However, instead of solving the non-trivial equation for surface pressure, in the model an approximate approach is used in which surface pressure is simply diagnosed from integrating density above a reference depth of 1500 m.

Additional tracers, including CFCs, dye and age tracers, have been added to the model as useful diagnostics. The air–sea fluxes of CFCs are computed following the Ocean Model Intercomparison Project (OMIP) protocol (Orr et al.2017).

Additional details of the ocean model are provided in Appendix B.

Figure 19(a) Historical global mean near-surface air temperature simulated by CLIMBER-X (black) compared to observations (blue) (Morice et al.2012) and CMIP6 models (grey). (b–d) Global mean near-surface air temperature in CLIMBER-X and CMIP6 models for idealized historical simulations with greenhouse gas concentration forcing only, natural (solar and volcanic) forcing only, and aerosol forcing only, respectively.

2.3 Sea ice model: SISIM

SISIM is a dynamic–thermodynamic sea ice model, with a single ice layer and a snow layer on top. SISIM also serves as a coupler between the atmosphere and ocean, and hence all surface energy fluxes over both sea ice and open ocean are computed in the sea ice model.

Figure 20Historical ocean heat content anomalies in (a) the top 700 m and (b) the top 2000 m simulated by CLIMBER-X (black) and derived from observations (blue) (Levitus et al.2012).

Sea ice thermodynamics are based on the zero-layer model of Semtner (1976). This involves the determination of accumulation and melting of snow/sea ice from above and accretion and melting of sea ice from below.

The surface energy balance equation is solved separately for the ice-free and ice-covered areas of each oceanic grid cell. Over sea ice, it is solved implicitly for the skin temperature, which is the temperature that balances all energy fluxes at the surface. The fluxes to be balanced include the net shortwave and longwave radiation at the surface, the sensible and latent heat fluxes to the atmosphere and the conductive heat flux into the snow/ice. Whenever the skin temperature over the snow/sea ice layer is above the melting point, the skin temperature is reset to 0C and the excess energy is used to melt snow/sea ice. Snowfall increases snow thickness, and if the snow depth exceeds 1 m, the snow excess is converted to sea ice. The net energy at the base of the ice layer is determined by the conductive heat flux through the snow/ice layer and the turbulent heat flux between the ice and the seawater below. Whether accretion or ablation occurs depends on the sign of the net energy at the ice–water interface. Over seawater, the heat flux into the ocean is computed as the residual of the radiation and surface energy fluxes, whereby the surface energy fluxes are computed using sea surface temperature. Sea ice forms whenever the top-layer ocean temperature drops below the freezing point of seawater.

Figure 21Equilibrium climate sensitivity versus transient climate response for the CLIMBER-X and CMIP6 models. CMIP6 model data are from Nijsse et al. (2020).

The surface sensible and latent heat fluxes are proportional to the temperature and moisture gradients between the surface and near-surface air, with the exchange coefficients depending on atmospheric stability through a bulk Richardson number. Sea ice albedo is temperature-dependent, with a decrease in albedo as the sea ice skin temperature approaches the melting point. This represents the decrease in albedo following the formation of meltwater ponds on the ice surface. The snow albedo scheme is the same as used in the land model and includes a dependence on snow grain size and dust and soot concentration following Dang et al. (2015). The conductive heat flux within the sea ice/snow layer is assumed to be directly proportional to the temperature difference between the surface skin and the salinity-dependent freezing temperature at the base of the ice layer (Semtner1976). The sea ice heat conductivity is scaled to represent sub-grid thickness distribution following Fichefet and Maqueda (1997). The turbulent heat flux at the ice base is determined from the difference between salinity-dependent freezing point temperature and top-layer ocean temperature using a constant exchange coefficient following McPhee (1992) and Weaver et al. (2001).

Figure 22Global feedback parameters for the CLIMBER-X and CMIP5 models. The feedbacks are (from left to right) Planck, water vapour, cloud, lapse rate, albedo, water vapour + lapse rate, temperature (Planck + lapse rate), and the sum of all feedbacks. CMIP5 data are from IPCC AR5.


Sea ice drift velocities are computed from the momentum balance equation (e.g. Hibler1979) with the elastic–viscous–plastic rheology of Hunke and Dukowicz (1997) and Bouillon et al. (2009). Numerically the solution of the momentum equation follows the implementation in GFDL model SIS2 (Adcroft et al.2019; Delworth et al.2006), discretized on the Arakawa C grid to match the ocean model grid. The derived velocities are then used to advect the grid cell mean ice and snow thicknesses and the sea ice concentration using a flux-corrected transport scheme (Zalesak1979). No explicit diffusion is applied.

Sea ice is allowed to cover the ocean fraction of the grid cell that is not occupied by shelf ice. That is, floating shelf ice restricts the domain available to sea ice.

Additional details of the sea ice model are provided in Appendix C.

Figure 23Zonal mean feedback parameters for CLIMBER-X. The total feedbacks (black solid lines) are further separated into contributions from longwave (LW, dashed black lines) and shortwave (SW, dotted black lines) radiation. Cloud feedbacks are additionally explicitly decomposed into LW and SW in panels (c) and (d), with the different colours representing feedbacks from changes in cloud fraction (blue), cloud optical thickness (red), and cloud top height (green).


2.4 Land model: PALADYN

CLIMBER-X includes the land model PALADYN (Willeit and Ganopolski2016), which serves as a land surface scheme for the exchange of energy and water between the surface and the atmosphere but also represents vegetation dynamics and, more generally, the land carbon cycle processes. The model treats in a consistent manner the interaction between atmosphere, terrestrial vegetation and soil through the fluxes of energy, water and carbon. Permafrost is treated explicitly, both in physical processes and as an important carbon pool. PALADYN distinguishes eight surface types, namely five plant functional types, bare soil, land ice and lakes. Over each surface type the model solves the surface energy balance and computes the fluxes of sensible, latent and ground heat and upward shortwave and longwave radiation. It includes a single snow layer. Vegetation and bare soil share a single soil column. The soil is vertically discretized into five layers where prognostic equations for temperature, water and carbon are consistently solved. Phase changes of water in the soil are explicitly considered. A surface hydrology module computes precipitation interception by vegetation, surface runoff and soil infiltration. The soil water equation is based on Darcy's law. Given soil water content, the wetland fraction is computed based on a topographic index. Modifications to the model relative to the description in Willeit and Ganopolski (2016) are described next.

Figure 24State dependence of global feedback parameters for CLIMBER-X. The feedbacks are for CO2 doubling starting from different initial CO2 concentrations (140, 280 and 560 ppm, from left to right). The feedbacks are computed from simulations where vegetation is prescribed at its pre-industrial state. The vegetation feedback is diagnosed from the difference in the sum of all feedbacks between simulations with dynamic and fixed vegetation. The sum of all the feedbacks also includes the contribution from the Planck feedback, which is not shown.


PALADYN has been complemented with a parameterization of dust emissions following the CLIMBER-2 scheme as described in Bauer and Ganopolski (2010). Dust is emitted into the atmosphere from deserts and grasslands whenever the soil is dry enough and the wind speed exceeds a critical threshold.

The snow albedo scheme has been refined with the inclusion of the effect of dust and soot on snow albedo following Dang et al. (2015), and the snow grain size parameterization has been retuned using output of the regional climate model MARv3.6 simulations for Greenland (Fettweis et al.2017). The parameterization for sub-grid snow cover fraction has been updated with a dependence on the standard deviation of orography following Roesch et al. (2001). More details on the updated surface albedo scheme can be found in Appendix D.

Figure 25Relative (a) global water content change, (b) global precipitation change and (c) salinity pattern amplification versus global temperature change for the CLIMBER-X and CMIP5 models from the 1 % yr−1 CO2 increase experiment. Each circle represents 1 year of the simulations.


Figure 26Zonal mean temperature response for atmospheric CO2 concentrations of (a) 180 ppm and (b) 560 ppm relative to 280 ppm with (solid) and without (dashed) vegetation feedback.


Figure 27Simulated dominant plant functional types at equilibrium for (a) the pre-industrial control, (b) 180 ppm and (c) 560 ppm of atmospheric CO2.

CLIMBER-X, and therefore also PALADYN, does not resolve the diurnal cycle. However, diurnal variations can be important, particularly for highly non-linear processes such as snowmelt. For instance, even if the daily mean surface air temperature is below freezing, the daily maximum temperature could still be substantially above the freezing point, thus leading to snowmelt during the day. In the presence of a relatively thin snow layer it is likely that the snow melted during the day would not completely refreeze during the night but would infiltrate into the soil or contribute to surface runoff. Following this line of thought, a parameterization of the diurnal cycle for snowmelt, partly following Krapp et al. (2017), has been implemented to accelerate the spring-time retreat of snow-covered area, which was too slow in the original PALADYN formulation (Willeit and Ganopolski2016). The diurnal cycle of skin temperature is assumed to have a sinusoidal shape with the amplitude depending on the difference between daily mean and daily minimum net shortwave radiation at the surface. The latter is computed from the diurnal cycle of insolation at the top of the atmosphere, which is available from the insolation module, using the daily mean surface albedo. If the resulting daily maximum skin temperature is above the freezing point, the energy flux available for snowmelt is computed following Krapp et al. (2017).

The possibility of prescribing land use changes has also been introduced into PALADYN following Burton et al. (2019). A fraction of each grid cell is prescribed as being used for agriculture, with no distinction between cropland and pasture being made. Land use is then represented as a limitation to the space available for a plant functional type (PFT) to expand into. For instance, the three woody PFTs (broadleaf trees, needleleaf trees and shrubs) are prevented from growing in the agricultural fraction, while the two grass PFTs (C3 grass and C4 grass) are allowed to grow anywhere in the grid cell and are simply interpreted as agricultural grasses if they grow into the pasture or cropland grid cell fraction. The representation of the effect of land use change on land carbon fluxes will be discussed in more detail in the CLIMBER-X companion paper about the global carbon cycle.

Since the original PALADYN publication, lakes have been introduced as an additional surface type in the model. Lake fractions are prescribed. Lake thermodynamics follows the implementation in CLM4.5 (Oleson et al.2013; Subin et al.2012) and includes freezing and melting of lake water and a single snow layer on top of lake ice. Temperatures and ice fractions are simulated for seven vertical layers, with a top-layer thickness of 1 m and the bottom-layer thickness depending on the depth of the lake, which is determined from lake bathymetry. For unfrozen lakes, convective mixing occurs in the lake if the density stratification becomes unstable. Additional vertical mixing occurs when partly frozen layers exist below not fully frozen layers in order to keep ice contiguous at the top of the lake. It is additionally always ensured that the lake is mixed to a minimum depth of 20 m. The surface energy fluxes over lakes are computed similarly to the other surface types as described in Willeit and Ganopolski (2016).

3 Model tuning

Tuning is an essential step in the development of any Earth system model. In such models a number of parameters are uncertain, either because they are poorly constrained by observations or because they are used to describe processes that cannot be explicitly represented in the model (e.g. Mauritsen et al.2012). These uncertain parameters are then adjusted in such a way as to minimize some measure of the distance between simulated and observed characteristics. In CLIMBER-X some processes that are explicitly represented in state-of-the-art climate models are parameterized instead, introducing some additional parameters that have to be properly calibrated. This refers in particular to the vertical structure of the atmosphere and friction in the ocean.

CLIMBER-X in general has been tuned to present-day climatological fields that are available from observations, reanalysis products and, in some cases, results from state-of-the-art model simulations. Climate feedbacks and model sensitivities have been considered an additional source of information in the process of model calibration, in particular for the parameterization of cloud properties.

The atmosphere model tuning involved a few steps. In the first step the longwave and shortwave radiative transfer schemes have been tuned offline. Shortwave and longwave radiation parameters were tuned in offline mode in order to minimize the root mean square error in top-of-the-atmosphere and surface radiation fluxes. Additional constraints on model parameters are provided by the total greenhouse effect separation presented in Schmidt et al. (2010) and by the radiative kernels of Shell et al. (2008), Block and Mauritsen (2013), Smith et al. (2018) and Pendergrass et al. (2018). More details on the offline tuning of shortwave and longwave radiation are available in Appendices A7 and A8. In the second step, the atmosphere model has been tuned in a coupled model set-up but using sea surface temperature relaxation for the ocean model in order to ensure that sea surface temperatures are close to present-day observations. Tuning parameters in the atmosphere model include parameters describing cloud cover fraction, cloud top height, cloud optical thickness, temperature lapse rate, height of the tropopause, strength of the Hadley and Ferrel cells and macro diffusivities for heat and water.

The land model has been tuned offline forced by prescribed climate variables as described in Willeit and Ganopolski (2016). The tuning of the snow grain size needed in the snow albedo parameterization has been performed offline using results of simulations with the regional climate model MAR for Greenland (Appendix D).

The ocean model tuning has been performed offline using sea surface temperature and sea surface salinity relaxation to observations. Ocean tuning parameters include friction and isopycnal and diapycnal diffusivities.

The sea ice model has been tuned in the fully coupled model set-up. Sea ice model tuning parameters include the thickness of newly formed sea ice, bare sea ice albedo and friction velocity below sea ice.

The coupling between the different model components introduces feedbacks that can substantially alter the performance of the model compared to the various offline or semi-coupled set-ups. A fine-tuning of the fully coupled model is therefore essential and has constituted a substantial part of the CLIMBER-X tuning procedure. The coupled model tuning was performed relying on a largely subjective and heuristic approach, with the goal of primarily minimizing the model–data mismatch in terms of global mean and zonal mean surface air temperature, global mean radiative fluxes at the top of the atmosphere (TOA) and at the surface and seasonal sea ice cover.

In the following we present results of model simulations performed using a single set of model parameters, produced based mainly on subjective ideas about “good” model performance. In the future we are planning to generate an ensemble of model versions selected using more formalized criteria to allow for a quantitative assessment of model uncertainties.

4 Model evaluation for the historical period and present day

The historical period, with its extensive availability of climate data from direct observations, forms the basis for the evaluation of any climate model. Here we evaluate the performance of CLIMBER-X for the climatological period from 1981 to 2010 and for the period of time from 1850 to 2015, corresponding to the historical period covered by CMIP (Coupled Model Intercomparison Project) simulations. For that we compare the results of CLIMBER-X simulations to different observation-based datasets as well as atmosphere and ocean reanalysis data. To give an overview of how CLIMBER-X compares to state-of-the-art general circulation models, we also include results from model simulations from the recent coupled model intercomparison projects. CMIP5 (Taylor et al.2012) and CMIP6 (Eyring et al.2016) model data are used interchangeably according to data availability.

The forcings for the historical CLIMBER-X simulations include variations in solar radiation (Matthes et al.2017), radiative forcing of volcanic eruptions (IPCC2013), globally uniform CO2, CH4 and N2O concentrations from Köhler et al. (2017), globally uniform CFC11 and CFC12 concentrations from Meinshausen et al. (2017), and 3-D O3 concentrations and 2-D SO4 load from the ensemble mean of CMIP6 models and land use change (pasture and cropland fractions) from Ma et al. (2020). The model is initialized from a 5000-year equilibrium simulation with pre-industrial boundary conditions.

4.1 1981–2010 climatology

In the following, different simulated climatological characteristics are compared to observations to assess the model performance for the present day. Unless stated otherwise, the comparison to observations is for the time interval from 1981 to 2010.

The global mean near-surface air temperature averaged over the time period 1981–2010 in CLIMBER-X is 14.1 C, which compares well to the 14.35 C in ERA-Interim reanalysis and is in the middle of the wide range of roughly 13–15 C spanned by CMIP5/6 models (e.g. Bock et al.2020). Global temperature is largely determined by the global radiation and energy budget, whose components are in good agreement with observations and CMIP5 models, as shown in detail in Table 1. The strength of the hydrological cycle is a bit overestimated in the model, resulting in higher evaporation and precipitation than observed, particularly over the ocean (Table 2). Note that, while the global latent heat flux in CLIMBER-X is underestimated compared to the observations in Table 1, the global evaporation is overestimated relative to the observational estimates in Table 2. This is a result of the use of different observation-based estimates for the energy fluxes and for the hydrological cycle and is thus ultimately a consequence of the uncertainty in these estimates.

Table 1Earth's radiation and energy budget in CLIMBER-X compared to observations and CMIP5 models from Wild et al. (2013). The unit is W m−2.

Download Print Version | Download XLSX

Table 2Global hydrological cycle in CLIMBER-X compared to observations (Trenberth et al.2007). The unit is 1015kg yr−1.

Download Print Version | Download XLSX

Atmosphere and ocean dynamics play an important role in the climate system by transferring energy from low to high latitudes, thereby reducing the Equator-to-pole temperature gradient arising from differential solar heating. The total meridional energy transport simulated by CLIMBER-X is compared to observation-based estimates in Fig. 2a. Overall the agreement is good, with a slight underestimate of the peak poleward energy transport. The meridional energy transport by the atmosphere, which constitutes the dominant contribution to total energy transport, matches well with observations (Fig. 2b). A separation of the different atmospheric processes contributing to latitudinal energy transport is also shown in Fig. 2b, including the partition between dry static energy and latent heat fluxes and the contribution by eddies. Eddies dominate the meridional energy transport at mid to high latitudes, while the Hadley cells are the most important contribution in the tropics, in agreement with a similar decomposition performed with other models (e.g. Yang et al.2015).

The mean meridional circulation in the atmosphere is characterized by the presence of six cells and is a defining feature of atmospheric dynamics. Figure 3 shows that the parameterization of the meridional circulation employed in CLIMBER-X does a reasonable job at describing the Hadley cells in the tropics and the Ferrel cells at mid latitudes.

The zonal mean SLP, which is directly related to the mean meridional circulation, shows good agreement with reanalysis data for both December–February (DJF) and June–August (JJA) (Fig. 4c, f). The model also reproduces the main features of observed azonal SLP, such as the highs over land and lows over the ocean in Northern Hemisphere winter (Fig. 4a, b) and vice versa during summer (Fig. 4d, e).

The simulated near-surface air temperature is compared to reanalysis in Fig. 5. In the zonal mean, the simulated temperatures fall mostly within the range given by different CMIP6 models (Fig. 5c, g). The main model biases include too warm temperatures over eastern continents in NH winter, too cold tropics, particularly over land, and a warm bias over East Antarctica (Fig. 5b, f). The absolute model bias is generally comparable to the bias seen in CMIP6 models, except for the tropics (Fig. 5d, h). The cold bias in the tropics is not only a surface feature, but is also persistent throughout the troposphere (Fig. 6). The rest of the 3-D temperature structure is well simulated by the model, indicating that the vertical variations in temperature can be reasonably well-described by the employed lapse rate parameterization. Zonal mean temperature biases are below a few degrees over large parts of the domain. During the winter months, CLIMBER-X also captures the near-surface temperature inversions at high latitudes. The simulated tropopause height shows a realistic latitudinal profile but is generally a few kilometres too low in the tropics (Fig. 6). CLIMBER-X also reproduces the higher seasonal variations in temperature over the continents compared to the ocean, in particular at high latitudes (Fig. 7).

The simulated atmospheric relative humidity is generally high in the planetary boundary layer and shows pronounced minima in the subtropics, broadly in agreement with observations (Fig. 8). The model does a reasonably good job at reproducing the observed precipitation distribution (Fig. 9). In terms of zonal mean precipitation, the peak associated with the intertropical convergence zone (ITCZ), the minima in the subtropics and the maxima at mid latitudes are well captured by the model (Fig. 9c, f). The main deficiencies are found in the subtropics, with too much precipitation simulated over the ocean in the subsidence areas (Fig. 9a, b, d, e). This is partly related to the too weak subtropical high-pressure systems (Fig. 4).

Cloud cover in CLIMBER-X shows the characteristic latitudinal profile, with minima in the subtropical subsidence areas and maxima in the tropics and at mid to high latitudes (Fig. 10a, f). A realistic simulation of cloud cover is a prerequisite for a good representation of radiation fluxes. Both shortwave and longwave radiation fluxes at the top of the atmosphere are in good agreement with satellite observations and reanalysis (Fig. 10b, d, g, i). The zonal mean radiative fluxes are generally within the range of CMIP5 models (Fig. 10c, e, h, j). Net shortwave radiation at TOA is slightly overestimated at high latitudes in NH summer (Fig. 10h), while net longwave TOA radiation exhibits some systematic biases in the tropics (Fig. 10e, j).

The ocean overturning circulation in CLIMBER-X is characterized by the presence of an Atlantic overturning cell and by Antarctic bottom water formation (Fig. 11). The maximum of the Atlantic overturning streamfunction at 26 N is 18.5 Sv, a bit higher than indicated by observations (Frajka-Williams et al.2021) (Fig. 12). The simulated Atlantic meridional circulation (AMOC) penetration depth of  3800 m, as measured by the zero crossing in the streamfunction, is about 500 m too shallow compared to that directly observed, a problem common also to many CMIP6 models (Fig. 12). The maximum meridional heat transport by the Atlantic Ocean is  1.12 PW, slightly lower than observation-based estimates of  1.25 PW (e.g. Johns et al.2011). The discrepancy between stronger-than-observed AMOC and the weaker-than-observed Atlantic meridional heat transport can be explained by biases in the simulated vertical structure of the transport (Fig. 12) and in the Atlantic Ocean temperature field (Fig.14). Deep water forms in the model at several locations in the northern North Atlantic, i.e. in the Labrador Sea south of Greenland and in the Nordic Seas, in agreement with observational estimates (Fig. 13). Deep water is also formed at several locations around Antarctica, mostly on the continental shelf, as shown by the annual maximum mixed-layer depth in Fig. 13. No deep water is formed in the North Pacific.

Ocean temperature and salinity fields compare well to observations in the deep ocean (Figs. 1415), with model biases mostly concentrated in the upper 1000 m. Biases in simulated temperature include too cold intermediate waters in the subtropics in the Atlantic and Indian oceans and too warm surface water in the North Pacific (Fig. 14c, f, i). Salinity biases of up to 1 psu are present in all ocean basins in the upper  1000 m (Fig. 15c, f, i).

The seasonality in sea ice area in both the NH and Southern Hemisphere (SH) is well reproduced by CLIMBER-X (Fig. 16) and is mostly within the range of CMIP6 models. The spatial extent of minimum and maximum sea ice cover is also in generally good agreement with observations (Fig. 17). Arctic winter sea ice cover is overestimated in the Fram Strait, while it is underestimated in the Sea of Okhotsk (Fig. 17a, e). Minimum and maximum sea ice extents in the Southern Ocean are well represented in the model (Fig. 17b, d, f, h).

The simulated total permafrost area is 18.5×106km2, close to the observed value of 18.8×106km2 (e.g. Tarnocai et al.2009). In terms of spatial extent, permafrost area is in good agreement with observations over Eurasia, while it is underestimated in eastern Canada (Fig. 18).

4.2 Simulations for the historical period

The CLIMBER-X-simulated historical evolution of global mean temperature is compared to observations (Morice et al.2012) and CMIP6 models in Fig. 19a. The model reproduces rather well the observed historical temperature trends and the response to volcanic eruptions. CLIMBER-X does not represent internal climate variability, and its results therefore cannot be compared one to one to observations. However, the simulated temperature shows a very good match with the ensemble mean of CMIP6 models, where internal variability has effectively been removed. The contribution of the different forcings to the historical temperature evolution also shows good agreement with the corresponding CMIP6 ensemble means, except for the last 2 decades, when CMIP6 models tend to overestimate the observed temperature change (Fig. 19b–d).

The rate of heat uptake by the ocean is consistent with observations (Levitus et al.2012) until around the year 2000 but is overestimated after that (Fig. 20). However, the historical ocean heat uptake in CLIMBER-X is in better agreement with observations compared to most EMICs, including CLIMBER-2 (Eby et al.2013).

5 Model sensitivities

Present-day observations provide a relatively poor constraint on model sensitivities to different climate forcings. Comparison to state-of-the-art general circulation models for experiments with different forcings and boundary conditions is therefore crucial for a model like CLIMBER-X. To this end we performed a comprehensive analysis of climate feedbacks, compared the response to changes in CO2 for standard CMIP abrupt4xCO2 ​​​​​​​ and 1 % yr−1 CO2 increase experiments, evaluated the vegetation feedback and tested the response to last glacial maximum boundary conditions, which provides insights into the model response to different orbital configurations, topographies and land–sea masks. We also performed standard freshwater hosing experiments to investigate the stability properties of the Atlantic meridional ocean circulation. An overview of the results of these experiments is presented next.

5.1 Transient climate response, climate sensitivity and Charney feedbacks

The equilibrium climate sensitivity of the standard version of CLIMBER-X is 3.3 K as computed from the temperature change at equilibrium (5000-year-long experiment) for a doubling of atmospheric CO2. This is in the middle of the range of 1.5–4.5 initially derived by Charney et al. (1979), which is also the estimated range given by the latest IPCC report (IPCC2021). The transient climate response of CLIMBER-X, defined as the global temperature change at the time of CO2 doubling in 1 % yr−1 CO2 increase experiments, is 1.9 K. A plot of equilibrium climate sensitivity versus transient climate response shows that CLIMBER-X is also well within the range of CMIP6 models (Fig. 21).

CLIMBER-X includes code to diagnose the strength of the different climate feedbacks, which allows for a more detailed analysis of the processes controlling climate sensitivity in the model. The feedbacks are evaluated using the partial radiation perturbation method (Bony et al.2006; Wetherald and Manabe1988). In this method, partial derivatives of model top-of-the-atmosphere radiation with respect to changes in modelled fields (such as water vapour, lapse rate and clouds) are determined by diagnostically rerunning the model radiation code.

The global feedback parameters for CLIMBER-X computed for CO2 doubling relative to 280 ppm are shown in Fig. 22 and generally compare well to feedback parameters computed for different models (e.g. Bony et al.2006). The water vapour feedback is the largest positive feedback in the model, followed by cloud and albedo feedbacks. The lapse rate feedback is globally negative, in agreement with CMIP5 models (Fig. 22).

A look at the zonal mean feedback parameters gives further insight into the spatial distribution of the feedbacks (Fig. 23). The albedo feedback is large at high latitudes and is related to the sea ice retreat and reduced snow cover in a warmer climate (Fig. 23e). The water vapour feedback is associated with an increase in water vapour content in a warmer atmosphere and is larger in the tropics than at the poles (Fig. 23f). The lapse rate feedback is negative in the tropics because of the larger warming of the mid–upper troposphere relative to the surface, while it is large and positive at high latitudes as a result of a pronounced erosion of surface temperature inversions mainly due to retreating sea ice (Fig. 23g). These results are all in agreement with feedbacks diagnosed in different general circulation models (e.g. Colman et al.2001; Crook et al.2011). Feedbacks related to clouds are the most uncertain and account for a large portion of the spread in climate sensitivity within current general circulation models (e.g. Zelinka et al.2020). Clouds affect both longwave and shortwave radiation at the top of the atmosphere through different processes. Cloud feedbacks in CLIMBER-X are shown in Fig. 23b–d, including a separation into contributions from changes in cloud fraction, cloud optical thickness and cloud height. The net cloud feedback is positive at all latitudes, with the notable exception of the Southern Ocean, where a pronounced increase in optical thickness with warming causes a large negative shortwave cloud feedback. Shortwave and longwave radiation cloud feedbacks generally act to (at least partly) compensate in CLIMBER-X (Fig. 23c–d), in agreement with ESMs (e.g. Zelinka et al.2012). The effects of cloud optical depth changes are generally larger in the shortwave, while changes in cloud height have a larger effect on longwave radiation. Combined, all feedbacks are mostly negative, except at high latitudes, where sea ice melting leads to large positive albedo and lapse rate feedbacks (Fig. 23h).

Climate feedbacks can be state-dependent. This is explored in CLIMBER-X by repeating the feedback analysis for a CO2 doubling but starting from different baseline CO2 concentrations, i.e. 140 and 560 ppm, in addition to the standard feedback analysis shown above, which was performed starting from a pre-industrial CO2 of 280 ppm. The albedo feedback shows a pronounced state dependence in CLIMBER-X, as shown in Fig. 24. In particular, the albedo feedback strongly increases in colder climates due to extended snow and sea ice cover. This is in agreement with previous studies, e.g. Colman and McAvaney (2009). The cloud feedback increases slightly with global temperature, while the sum of water vapour and lapse rate feedbacks decreases as climate warms.

Hydrological sensitivity, which quantifies the relative global precipitation change per unit change in global temperature, is an important measure of the response of the hydrological cycle to climate change (e.g. Held and Soden2006). In agreement with CMIP models, CLIMBER-X shows an increase in global precipitation by  2 % per degree global temperature increase in the 1 %-per-year CO2 increase experiment (Fig. 25b), while the atmospheric water content increase approximately follows the Clausius–Clapeyron increase of 7 % K−1 (Fig. 25a). It has also been observed that the sea surface salinity pattern is amplified under global warming, with salinity increasing in high-salinity regions and decreasing in low-salinity regions due to an intensification of existing patterns of precipitation–evaporation (e.g. Durack et al.2012). The salinity pattern amplification (as defined by e.g. Durack et al.2012) in the 1 %-per-year CO2 increase simulation in CLIMBER-X is  5 % per degree global warming (Fig. 25c), in good agreement with estimates from Zika et al. (2018) for the historical period. The salinity pattern amplification is therefore larger than the hydrological sensitivity of the model, possibly due to the effect of ocean warming on surface salinity (Zika et al.2018).

5.2 Vegetation feedback

Changes in vegetation structure and its spatial distribution can affect the climate through the effect on surface energy and water fluxes (e.g. Levis et al.1999; Bala et al.2006; Falloon et al.2012). In CLIMBER-X, the vegetation feedback amplifies temperature change by more than 1 C at mid to high northern latitudes for a reduction of CO2 to 180 ppm (Fig. 26a), while it is small for CO2 doubling (Fig. 26b). A similar asymmetry in the vegetation feedback between low and high CO2 has also been found in CLIMBER-2 (Willeit et al.2014). The reason for the strong positive vegetation feedback for low CO2 originates from a pronounced southward retreat of boreal forest (Fig. 27b), which causes a substantial increase in surface albedo through the missing snow-masking effect of trees (e.g. Bonan2008). In CLIMBER-X, a CO2 increase leads in general to an expansion of forests and a reduction in desert area, while the opposite happens for lower CO2 levels (Fig. 27). The strong state dependence of the vegetation feedback can also be seen in Fig. 24, with a very strong feedback for a CO2 concentration of 140 ppm and negligible vegetation feedback for climates warmer than the present day.

5.3 Last glacial maximum

For the last glacial maximum (LGM), we prescribe boundary conditions following the PMIP4 protocol (Kageyama et al.2017), with the GLAC-1D ice sheet, bathymetry and land–sea mask reconstruction (Tarasov et al.2012). The simulation is started from the present-day equilibrium followed by a switch to LGM boundary conditions during which the topography, bathymetry and ocean volume are adjusted. Total ocean salinity is conserved in this process. The model is then run for 5000 years to ensure equilibrium is reached.

The simulated global cooling at the LGM relative to the pre-industrial is 6.2 K, similar to the cooling produced by CLIMBER-2 (Ganopolski et al.1998). This matches the most recent reconstruction-based estimate by Tierney et al. (2020); however, it is on the cold side of range produced by PMIP4 models (3.3–7.2 K) (Kageyama et al.2021). The vegetation feedback in the model is responsible for  0.5 K of the simulated cooling. The zonal mean annual temperature change is compared to PMIP3/4 models in Fig. 28. CLIMBER-X shows a cooling in the tropics that is more pronounced than in most other models, while at high latitudes the temperature difference falls well inside the PMIP3/4 range. In terms of sea surface temperatures, the model results agree well with the proxy-based reconstruction by Tierney et al. (2020) (Fig. 29), which also shows a pronounced cooling in the tropics as opposed to e.g. Paul et al. (2021) (Fig. 29c).

Figure 28Last glacial maximum annual mean zonal near-surface air temperature differences relative to pre-industrial compared to PMIP3/4 models (Kageyama et al.2021).

Figure 29Last glacial maximum annual mean sea surface temperature differences relative to pre-industrial compared to reconstructions (Paul et al.2021; Tierney et al.2020).

The global ocean cools by 2.4 C, in good agreement with 2.57 ± 0.24 C in Bereiter et al. (2018), with the deep ocean temperature being close to the freezing point of seawater (<-1C) in all ocean basins, in accordance with Adkins et al. (2002). The AMOC is weaker and shallower at the LGM relative to the pre-industrial, with a maximum overturning strength of  14 Sv at 26 N and extending down to a depth of  2500 m (Fig. 30a). This is contrary to most PMIP3/4 models, which tend to produce a more vigorous and deeper Atlantic overturning at the LGM (Weber et al.2007; Kageyama et al.2021) but is possibly in better agreement with proxy reconstructions (e.g. McManus et al.2004; Bohm et al.2015). The AMOC weakening is accompanied by a strengthening of Antarctic bottom water formation (Fig. 30b), which is related to an increase in brine rejection associated with a pronounced expansion of sea ice in the Southern Ocean (Fig. 31b), similarly to what has been found in other models (e.g. Nadeau et al.2019; Shin et al.2003; Stouffer and Manabe2003). The large simulated sea ice expansion in the Southern Ocean is also largely consistent with the latest proxy reconstructions by Lhardy et al. (2021) (Fig. 31b).

Figure 30Ocean overturning circulation at the last glacial maximum: (a) global and (b) Atlantic.


Figure 31Difference in seasonal sea ice area in the (a) NH and (b) SH between the last glacial maximum and the pre-industrial. CLIMBER-X results (black) are compared to CMIP5 models (grey). In panel (b) the seasonal minimum and maximum changes in SH sea ice area from proxy estimates of Lhardy et al. (2021) are also indicated.

5.4 AMOC stability

The meridional overturning circulation in the Atlantic Ocean plays an important role in the global climate system. Since the pioneering work of Stommel (1961), who used a simple box model to show that the AMOC could have multiple stable states, the stability of the AMOC has received considerable attention, with several studies using ocean general circulation models confirming the bi-stable nature of the system (e.g. Manabe and Stouffer1988; Rahmstorf1995).

Consistent with findings by other EMICs (Rahmstorf et al.2005) and low-resolution atmosphere–ocean GCMs (Hawkins et al.2011), CLIMBER-X also shows a hysteresis behaviour of the AMOC when the freshwater balance of the North Atlantic is perturbed (Fig. 32). The width of the hysteresis for a rate of change of 0.05 Sv kyr−1, which was also used by Rahmstorf et al. (2005), is  0.25 Sv. However, when quasi-equilibrium simulations (5000 years long) are performed with the model starting from AMOC on and off states for a given constant freshwater hosing flux, the width of the hysteresis is reduced by about a factor of 2 (Fig. 32). This indicates that a rate of change of 0.05 Sv kyr−1 is too high to allow a precise tracking of the actual AMOC stability diagram. Under pre-industrial conditions, the AMOC is in a mono-stable regime in the model although relatively close to the bi-stable regime. The present-day stability of the AMOC is still debated (e.g. Weijer et al.2019).

Figure 32AMOC hysteresis for freshwater hosing applied to the latitudinal belt between 20 and 50  N in the Atlantic. The hysteresis curve for experiments where the freshwater hosing is changed at a rate of 0.05 Sv per 1000 years is shown by the solid lines. The symbols and the connecting dotted lines indicate the AMOC strength for quasi-equilibrium simulations with a prescribed constant freshwater hosing rate initialized from an AMOC on state (red circles) and from an AMOC off state (blue crosses).


The different climatic conditions associated with the two equilibrium states of the AMOC are illustrated by mean annual surface air temperature and precipitation differences in Fig. 33. The temperature difference shows the classic seesaw pattern, with cooling in the NH and warming in the SH as a response to AMOC shutdown (Fig. 33a). The cooling reaches up to 10 K in the North Atlantic and is compensated by a warming of up to 10 K in the Southern Ocean. The cooling produced in the North Atlantic is consistent with GCM simulations (e.g. Vellinga and Wood2002; Jackson et al.2015; Pedro et al.2018), while the large warming in the Southern Ocean is not seen in these GCMs, possibly because they are not run into equilibrium. Precipitation changes in the AMOC off state compared to the on state are clearly seen in the tropics as a result of a pronounced southward shift in the intertropical convergence zone (Fig. 33b). Precipitation is also reduced in the North Atlantic, simply as a consequence of the colder climatic conditions.

Figure 33Annual mean (a) near-surface air temperature and (b) precipitation differences between AMOC off and on states.

6 Applicability and limitations of CLIMBER-X

CLIMBER-X does not resolve synoptic variability and does not exhibit interannual internal variability and is therefore not suited to investigating weather extremes or internal climate oscillations like ENSO. The atmospheric component of CLIMBER-X is based on a statistical–dynamical approach, which employs a number of significant simplifications and assumptions. Such an approach allows us to develop a model that is several orders of magnitude faster than GCMs with similar resolution; however, these simplifications and a set of parameterizations explicitly derived from present-day climate limit the model's applicability to climate states fundamentally different from the present.

7 Conclusions

We have described the major features of the climate component of the newly developed CLIMBER-X Earth system model. CLIMBER-X relies on the geostrophic approximation for the description of both atmospheric and oceanic circulation. It does not explicitly resolve atmospheric and oceanic eddies, but the effect of both is parameterized as a diffusive process. The simplified dynamics, together with the use of a daily time step for most processes and the lower spatial resolution of 5×5, provides a substantial advantage in terms of computational costs relative to general circulation models, which are based on the primitive equations. Nevertheless, in terms of the number of physical processes which are represented in the model, CLIMBER-X is comparable to state-of-the-art climate models. This is highlighted for instance by the realistic representation of climate feedbacks in the model. Moreover, CLIMBER-X also includes a model for the global carbon cycle and is coupled to an ice sheet model (both of which will be described in detail in forthcoming papers) and thus allows simulations of the evolution of the Earth system as a whole and the investigation of complex interactions and feedbacks between different components of the Earth system on timescales ranging from decades to hundreds of thousands of years. CLIMBER-X is therefore an ideal tool to explore the long-term past and future evolution of the Earth system.

Appendix A: SESAM detailed description

A1 Vertical structure

SESAM is based on a number of simplifications compared to current state-of-the-art atmospheric models. One such simplification is that, in all equations apart from that describing atmospheric dynamics, the atmospheric pressure at each model level is not computed using the hydrostatic approximation but is assumed to be an exponential function of the height above sea level only:

(A1) p ( z ) = p 0 e - z H a ,

where Ha=RdT0/g is the atmospheric height scale. The mean sea-level pressure p0 is determined from the atmospheric mass conservation condition:

(A2) Ω p ( z s ) d ω = M a g A e ,

where Ω indicates the surface of the Earth, zs is surface elevation above sea level, Ma is the total atmospheric mass and Ae is the area of the Earth's surface. Hence, the mean sea-level pressure will change with changing topography.

For atmospheric dynamics the dependence of pressure on atmospheric temperature is explicitly accounted for through the parameterization of sea-level pressure and through the thermal wind equation (see Appendix A2 below).

Air density is described as a function of elevation:

(A3) ρ ( z ) = ρ 0 e - z H a ,

with the reference density at sea level computed from the ideal gas law:

(A4) ρ 0 = p 0 R d T 0 .

The vertical profile of temperature is computed as follows:

(A5) T ( z ) = T a - z s z Γ ( z ) d z ,

where Ta is the prognostic atmospheric temperature at the surface and Γ is the lapse rate. Although the atmospheric temperature lapse rate is remarkably close to a constant value of  6.5 K km−1, the deviations from this value are important in order to reproduce both a realistic present-day climate and climate feedbacks. In SESAM, Γ is parameterized as

(A6) Γ ( z ) = Γ s , z z s + H Γ , s , Γ b + Γ t - Γ b z H Γ , t , z > z s + H Γ , s and z H T , 0 , z > H T .

In general, the lapse rate therefore linearly increases with height, with Γb and Γt depending on atmospheric humidity qa only:


In a layer close to the surface, the lapse rate depends on near-surface stability:

(A9) Γ s = c 4 Γ max 0 , T a - T , ocean , c 5 Γ T a - T , land and T a - T > 0 , c 6 Γ T a - T , land and T a - T < 0 , c 5 Γ T a - T , ice ,

where T is the skin temperature and Γs is limited to be lower than 7.5×10-3K m−1 over ocean and 10×10-3K m−1 over land and ice. In particular, Eq. (A9) allows SESAM to reproduce near-surface inversions which are important for surface climate.

The tropopause height HT is derived assuming that the stratosphere is in radiative equilibrium:

(A10) H T t = - c 1 tp R str , net + S .

The net radiation in the stratosphere, Rstr,net, includes the balance of longwave radiation and the shortwave radiation absorbed by ozone. The effect of atmospheric dynamics on tropopause height is explicitly included in the prescribed latitudinal profile of S, which only depends on the position (ϕITCZ) and width (ΔϕHad) of the Hadley cells:

(A11) S = c 2 tp 1 - c 3 tp 1 - sin 8 0.85 ϕ - ϕ ITCZ 0.5 Δ ϕ Had .

Potential temperature, which is a conserved quantity for adiabatic motions, is computed as

(A12) θ ( z ) = T ( z ) + Γ d z ,

where Γd=g/cp is the dry adiabatic lapse rate and cp is the specific heat capacity of air at constant pressure.

In CLIMBER-2 specific humidity was specified as decaying exponentially with height, but this can imply unrealistic relative humidities in the upper troposphere, where humidity is important for longwave radiation. This problem is also highlighted by the fact that there is no tight coupling between water vapour and lapse rate feedbacks in CLIMBER-2, contrary to what is observed in most climate models. To overcome this limitation, in SESAM the vertical profile of specific humidity is expressed through temperature and a parameterization of the relative humidity variation with height:

(A13) r ( z ) = r a z z pbl , r a e - z - z pbl H r z > z pbl and z z s + c 4 r , r a e - z s + c 4 r - z pbl H r z > z s + c 4 r and z H T , r st z > H T ,

where zpbl=zs+c5r is the elevation of the planetary boundary layer. The relative humidity height scale Hr is constant and uniform in the extratropics and depends on vertical velocity at 700 hPa in the tropics:

(A14) H r = f trop c 1 r e c 2 r w 700 + ( 1 - f trop ) c 1 r c 3 r ,

with ftrop=1-sin8φ, where φ is defined below in Appendix A2. Specific humidity is then computed as

(A15) q ( z ) = r ( z ) q sat ( T ( z ) , p ( z ) ) ,

where the specific humidity at saturation qsat is computed assuming saturation over ice at temperatures below −15C, saturation over water at temperatures above T0 and a weighted mean of saturation over water and ice in the intermediate temperature range.

Table A1Parameters for vertical structure.

Download Print Version | Download XLSX

A2 Dynamics

The dynamics of the atmosphere in SESAM is similar to that in CLIMBER-2 but with several notable improvements.

Horizontal velocity in the atmosphere is computed as the sum of geostrophic and ageostrophic components:

(A16) u = u g + u a .

The geostrophic components of velocity at any height within the troposphere are obtained using the thermal wind approximation:


where psl is the sea-level pressure, f is the Coriolis parameter, Re is the radius of the Earth and λ and ϕ are longitude and latitude, respectively. The ageostrophic wind components in the planetary boundary layer are computed from sea-level pressure and cross-isobar angle, α, as


As in Petoukhov et al. (2000) the ageostrophic wind in the planetary boundary layer is compensated in the upper troposphere in order to conserve mass in the atmospheric column. Since the geostrophic approximation is not valid close to the Equator, the Coriolis parameter is limited to |f|>3×10-5s−1 in the geostrophic and |f|>1×10-5s−1 in the ageostrophic wind equations.

The cross-isobar angle is determined from the condition that the shear stress is continuous between the Ekman layer and the surface layer (Petoukhov et al.2000):

(A21) C D U s 2 = U s ϵ sin α 2 | f | K v ,

where ϵ=1-sin2α, Us is the module of the surface wind, Kv is the kinematic vertical viscosity coefficient in the planetary boundary layer and the drag coefficient CD is

(A22) C D = κ ln z ref z 0 + z oro 2 ,

where κ is the von Karman constant, z0=100m a reference height, z0 the surface roughness length and zoro the orographic roughness, computed from the sub-grid-scale standard deviation of orography as

(A23) z oro = 0.004 σ oro .

Equation (A21) is solved for α with the approximation Us2K.

The near-surface wind components are computed using the Taylor model (see e.g. Hansen et al.1983) with the addition of a simple representation of katabatic winds:


Katabatic winds (uk, vk) are important in surface inversion conditions over slopes, as is at present mainly the case over the large Greenland and Antarctic ice sheets. They are included in the model based on a simple balance of buoyancy force and friction, ignoring Coriolis and background pressure gradients, following the Prandtl 1942 model (e.g. Fedorovich and Shapiro2009):


with h=100m and T2 m the near-surface air temperature.

Sea-level pressure is computed as the sum of zonally averaged and azonal components:

(A28) p sl = p sl + p sl * .

The zonal component of sea-level pressure is computed using a parameterization similar to that described in Petoukhov et al. (2000). This parameterization is based on the assumption that the strength and position of the major cells of atmospheric meridional overturning circulation are controlled by average meridional temperature gradients, the zonally averaged surface drag and surface elevation. The zonally averaged sea-level pressure is defined through the zonal mean meridional (ageostrophic) wind component in the PBL:

(A29) p sl ϕ = - v a ( ϕ ) ρ 0 f R e sin α cos α ,

where the zonally averaged meridional (sea-level) ageostrophic velocity is parameterized as

(A30) v a ( ϕ ) = - C i Δ T i j F z ( ϕ ) sin φ , ( i - 1 ) π < φ < i π ,


(A31) φ = 6 D had ϕ - ϕ ITCZ c 1 mmc ϕ - ϕ ITCZ 2 + 1 .

Ci are empirical parameters, i=1 corresponds to the Hadley cell, i=2 to the Ferrel cell and i=3 to the polar cell, and j=1 corresponds to the Northern Hemisphere and j=2 to the Southern Hemisphere. In contrast to CLIMBER-2, the same Ci values are used for the Northern Hemisphere and Southern Hemisphere. The position of the intertropical convergence zone, ϕITCZ, depends on the temperature difference between the two hemispheres:

(A32) ϕ ITCZ = c 2 mmc T NH - T SH .

Dhad controls the width of the Hadley cells and is a function of tropical temperature:

(A33) D had = c 3 mmc T trp - c 4 mmc .

This ensures that the Hadley cells are expanding with warming, consistent with empirical evidence (e.g. Hu et al.2018) and models (Frierson et al.2007). Tij is the mean temperature of each cell. The temperature gradients are proportional to meridional differences in zonal mean sea-level temperature:

(A34) Δ T i j = T ( ϕ i j ) - max ( T ( ϕ ) ) , i = 1 , T ( ϕ i j ) - T ( ϕ i - 1 j ) , i = 2 , 3 .

The latitudes used to compute the gradients are fixed at |ϕ1j|=π/6, |ϕ2j|=π/3 and |ϕ3j|=π/2. The topography factor is defined as

(A35) F z ( ϕ ) = 1 - z s / c 5 mmc .

The azonal sea-level pressure is the sum of thermal and orographic components:

(A36) p sl * = p sl , T * + p sl , O * .

The azonal sea-level pressure component arising from the thermally induced stationary planetary waves is computed following Petoukhov et al. (2000) and Petoukhov et al. (2003):

(A37) p sl , T * = - g p 0 H 0 2 R d T 0 2 T sl * ,

where Tsl is the skin temperature reduced to sea level using a constant lapse rate of 6.5 K km−1. The effect of topographic stationary waves is accounted for using the simple 1-D barotropic model for forced topographic Rossby waves of Charney and Eliassen (1949) as described in Held (1983) and Holton (2004), with the linearized vorticity equation written as

(A38) u ζ λ + β v + ζ τ e = - f H T 0.4 u z s λ ,

where ζ is the relative vorticity, β is the meridional derivative of the Coriolis parameter f, τe is the damping timescale due to Ekman pumping and u is taken to be the zonal mean zonal wind at 500 hPa. A meridional wave number corresponding to a latitudinal half-wavelength of 35  is assumed. The equation is solved independently for each latitudinal belt by Fourier expansion of the topography zs(λ) by writing the geostrophic wind vector and the relative vorticity in terms of a streamfunction Ψ. The FFTW3 library (Frigo and Johnson2005) is used for the Fourier transform. The sea-level pressure perturbation due to topographic stationary waves is then derived as

(A39) p sl , O * = Ψ f ρ ( 500 hPa ) .

Table A2Parameters for atmospheric dynamics.

Download Print Version | Download XLSX

A3 Thermodynamics

The energy balance equation, vertically integrated from the surface to the top of the atmosphere, is written as

(A40) Q T t = - 1 R e cos ϕ × λ z s H T ρ ( u θ + u θ ^ ) d z + ϕ z s H T cos ϕ ρ ( v θ + v θ ^ ) d z + c v - 1 SW a + LW a + L e P w + L s P s + SH ,

where QT=zsHTOAρTdz is the heat content of the atmospheric column, SWa is the shortwave radiation absorbed by the atmosphere, LWa is the net atmosphere longwave radiation balance, Pw is rainfall and Ps is snowfall, Le is the latent heat of evaporation and Ls the latent heat of sublimation, SH is the surface sensible heat flux, and cv is the heat capacity of air at constant volume. The horizontal heat transport due to synoptic processes, uθ^ and vθ^, is represented as macroturbulent diffusion as described in Appendix A5. In CLIMBER-2 only the non-thermal wind is used in the energy balance equation, with the beta effect accounted for separately. This allowed a large time step of up to 1 d. Because the CLIMBER-X horizontal resolution is much higher than in CLIMBER-2, a relatively short time step of  2 h is needed, and with such a small time step the energy equation is stable even using the full wind vector. The energy balance equation is solved for Ta and the near-surface air temperature is diagnosed as

(A41) T 2 m = T a + T 2 .

T2 m is also the temperature that is used to compute the surface turbulent sensible heat flux.

A4 Hydrology

The water balance equation, vertically integrated from the surface to the top of the atmosphere, is written as

(A42) Q q t = - 1 R e cos ϕ × λ z s H T ρ ( u q + u q ^ ) d z + ϕ z s H T cos ϕ ρ ( v q + v q ^ ) d z + E - P ,

where Qq=zsHTOAρqdz is the water vapour content of the atmospheric column, E is surface evaporation and P is total precipitation. The horizontal moisture transport due to synoptic processes, uq^ and vq^, is represented as macroturbulent diffusion as described in Appendix A5.

The water balance equation is solved for qa, and the near-surface air-specific humidity, which is used to compute the surface latent heat flux, is diagnosed as

(A43) q 2 m = r 2 m q sat ( T 2 m ) ,

with r2m=(ra+r)/2 and r=qa/qsat(T). Precipitation in the model is generated as follows:

(A44) P = max 0 , C + C slope + E r a r a max + Q q r a τ p ,

where C is moisture convergence into the atmospheric column by advection and diffusion and Cslope explicitly represents an additional moisture convergence due to synoptic activity on slopes. The first term on the right-hand side of Eq. (A44) represents a gradual removal of converging moisture by precipitation as atmospheric relative humidity increases, with all water entering an atmospheric column being added to precipitation when the atmospheric relative humidity reaches a maximum value ramax. The second term on the right-hand side of Eq. (A44) generates additional precipitation by removing atmospheric moisture with a given timescale τp and is applied only over land. The moisture convergence due to synoptic activity on slopes is computed assuming that moisture convergence is proportional to the vertical velocity induced by synoptic winds impacting on the slope:

(A45) C slope = c slope p K z s ρ 0 q a ,

where the synoptic wind is expressed through the eddy kinetic energy, K.

Table A3Parameters for the hydrological cycle.

Download Print Version | Download XLSX

A5 Synoptic processes

Horizontal fluxes of energy and water originating from unresolved synoptic variability are represented as a macroturbulent diffusion process (Petoukhov et al.2000):


The diffusivities AT and Aq are isotropic, are vertically uniform and depend on eddy kinetic energy K:


The different dependence of the diffusivities on eddy kinetic energy follows from Caballero and Hanley (2012).

The kinetic energy of synoptic eddies is described by an evolution equation:

(A52) K t = - u K + A T K + P K - D K .

Kinetic energy production is proportional to the maximum Eady model baroclinic growth rate (e.g. Hoskins and Valdes1990):

(A53) P K = c 1 syn + c 2 syn f N u z ,

where the Brunt–Vaisala frequency N is defined as

(A54) N = g θ θ z ,

and all vertical gradients are computed between 850 and 500 hPa. The dissipation of eddy kinetic energy is given by

(A55) D K = c 3 syn + c 4 syn C D K 3 / 2 .

The synoptic component of near-surface wind is computed from eddy kinetic energy as follows:

(A56) U syn = c 7 syn ϵ cos α K .

The synoptic vertical velocity at 700 hPa is calculated as

(A57) w syn = c 8 syn K .

The module of surface wind, which is used for calculation of the turbulent surface fluxes, is defined as

(A58) U s = u s 2 + v s 2 + U syn 2 ,

and the wind stress over the ocean is computed as


Table A4Parameters for synoptic processes.

Download Print Version | Download XLSX

A6 Clouds

Total cloud cover fraction in SESAM is a combination of clouds related to atmospheric relative humidity and vertical velocity (fcldr) and clouds related to surface temperature inversion conditions (fcldlow):

(A61) f cld = 1 - 1 - f cld r 1 - f cld low .

The relative humidity- and vertical-velocity-mediated cloud fraction, which provides the main contribution to cloud cover in the model, is given by

(A62) f cld r = c 1 cld + c 2 cld tanh c 3 cld w eff r a c 4 cld .

It is therefore proportional to raccld4, with the proportionality factor depending on the effective vertical velocity at cloud level (weff). In addition to the mean vertical velocity, weff also includes contributions from vertical velocities resulting from synoptic disturbances and sub-grid-scale orography:

(A63) w eff = w ( 700 hPa ) + c weff w syn + w oro .

The synoptic vertical velocity is given by Eq. (A57), while the orographic component is a function of surface wind speed and grid cell standard deviation of surface elevation (σoro):

(A64) w oro = c woro U s σ oro .

The fraction of low clouds related to the presence of surface temperature inversion, when r>ra, is defined as

(A65) f cld low = c 5 cld f freezedry ( r - r a ) + c 6 cld 2 c 6 cld r a c 4 cld .

The factor ffreezedry represents the freeze-dry mechanism following Vavrus and Waliser (2008) and decreases the low cloud amount in very cold and dry conditions,

(A66) f freezedry = 0.1 + 0.9 q a c 7 cld ,

and is limited to the range [0,1].

The cloud base height is assumed to coincide with the top of the planetary boundary layer, Hpbl. Cloud top height follows the height of the tropopause, modified by a factor depending on the mean vertical velocity at 700 hPa:

(A67) H cld = c 1 hcld + c 2 hcld H T 1 + c 3 hcld w ( 700 hPa ) .

Cloud optical thickness is parameterized as a function of surface air temperature, cloud fraction and column water content,

(A68) τ cld = c 3 τ 1 + tanh - T 2 m - T 0 - c 1 τ c 2 τ f cld Q q c 4 τ ,

and is further modified to account for the indirect effect of sulfate aerosols following Bauer et al. (2008).

Table A5Cloud parameters.

Download Print Version | Download XLSX

A7 Shortwave radiation

The computation of shortwave radiation fluxes is based on a two-stream delta-Eddington approximation of the transport equation in a gas–aerosol atmosphere (see Petoukhov et al.2003, for a more detailed description). The net downward shortwave radiation fluxes are computed at the top of the atmosphere and at the surface as a weighted mean of clear-sky (cs) and cloudy (cld) fluxes:


The clear-sky and cloudy fluxes are computed separately for two spectral bands – visible (vu) and near-infrared (ir) – and for each macro surface type. The individual components of the shortwave radiation flux are computed as


where fvu is the fraction of solar radiation in the visible and ultraviolet spectral range, αatm is the planetary albedo and Iatm is the atmosphere integral transmission function. Planetary albedos used to compute the net shortwave fluxes at TOA are defined as


The atmospheric scattering albedo depends on the cosine of the solar zenith angle (μ) and the optical thickness (τaer) and imaginary part of the refractive index (Raerim) of the aerosol load:


Cloud albedo is computed from cloud optical thickness as follows:


The atmospheric integral transmission functions used to derive the net shortwave fluxes at the surface are calculated as


The first terms on the right-hand side represent the direct effect on radiation reaching the surface, while the second terms represent the effect of multiple reflections between the surface and atmospheric scatterers and clouds.

The integral transmission functions for water vapour, aerosols, ozone and clouds are given by


where Mw and Maer are the effective absorber mass of water and aerosols in the column, respectively. In general, Mw and Maer are different for clear sky and cloudy sky and for the top of the atmosphere and the surface:


W is the column water content in g cm−2 and μ0 is the effective cosine of the zenith angle for diffuse radiation, Hq is the humidity height scale, Dcld is the geometrical thickness of clouds, fexp1=exp(-Hcld/Hq)-exp(-(Hcld+Dcld)/Hq) and fexp2=1-exp(-Hcld/Hq)/μ0.

Shortwave radiation parameters are tuned in offline mode using prescribed column water content, surface albedo and total cloud cover fraction from ERA-Interim reanalysis (Dee et al.2011) and cloud top height and cloud optical thickness from the ISCCP climatology (Rossow and Schiffer1999). The target is to minimize the root mean square error in top-of-the-atmosphere and surface shortwave radiation fluxes for clear sky and total sky. The resulting zonal mean net top-of-the-atmosphere shortwave radiation is shown in Figs. A1 and A2 for clear sky and total sky, respectively.

Additional constraints on model parameters are provided by the shortwave radiative kernels at the top of the atmosphere and at the surface, which are available for several climate models (Shell et al.2008; Block and Mauritsen2013; Smith et al.2018; Pendergrass et al.2018). The TOA kernels computed using the CLIMBER-X shortwave radiation code are compared to those available from the different climate models in Fig. A3. The sensitivity of TOA shortwave flux to changes in surface albedo is well captured by the shortwave radiation scheme of CLIMBER-X (Fig. A3a), while the sensitivity to water vapour changes is too large in the tropics (Fig. A3b). Note, however, that changes in water vapour have a generally small effect on TOA shortwave radiation, so that this bias is expected to have only minor implications for model results.

Table A6Shortwave radiation parameters.

Download Print Version | Download XLSX

Figure A1Zonal mean net shortwave radiation at the top of the atmosphere for clear-sky conditions as computed by the CLIMBER-X radiation code in the offline set-up compared to ERA-Interim reanalysis (Dee et al.2011) and CERES satellite observations (Loeb et al.2018). Panels (a) and (b) are for December–January–February and panels (c) and (d) are for June–July–August. Panels (a) and (c) show absolute values, while panels (b) and (d) show model biases.

Figure A2Zonal mean net shortwave radiation at the top of the atmosphere as computed by the CLIMBER-X radiation code in an offline set-up compared to ERA-Interim reanalysis (Dee et al.2011) and CERES satellite observations (Loeb et al.2018). Panels (a) and (b) are for December–January–February and panels (c) and (d) are for June–July–August. Panels (a) and (c) show absolute values, while panels (b) and (d) show model biases.

Figure A3Shortwave radiative kernels at the top of the atmosphere for (a) a 1 % increase in surface albedo and (b) water vapour increase corresponding to a 1 K warming of the atmosphere. The kernels computed using the CLIMBER-X radiation code are compared to kernels available from different climate models (Shell et al.2008; Block and Mauritsen2013; Smith et al.2018; Pendergrass et al.2018). Solid lines are for total-sky conditions, while dashed lines are for clear-sky conditions.

A8 Longwave radiation

Longwave radiative transfer follows in part from CLIMBER-2 as described in Petoukhov et al. (2003). Improvements relative to CLIMBER-2 include

  • relaxing the assumption that clouds act as black bodies for longwave radiation, which is not strictly true for thin clouds in cold regions,

  • a specific humidity profile derived from temperature and relative humidity, giving more consistent humidity/temperature values in the upper troposphere,

  • CH4, N2O, CFC11 and CFC12 being considered through the use of an equivalent CO2 concentration, and

  • extensive retuning in an offline set-up using observations and output of atmospheric GCMs, including radiative kernels.

The computation of longwave radiative transfer is based on the two-stream approximation, where the downward and upward longwave radiative fluxes at any height z are computed as


where B(z)=σT4(z) is the longwave emission at height z, Bs=ϵσTs4 the long wave emitted by the surface with an emissivity ϵ, σ the Stefan–Boltzmann constant, and D(z,z) the integral transmission function of the atmospheric layer confined between levels z and z. Longwave radiation fluxes are computed separately for clear sky and cloudy sky using the following integral transmission functions, which explicitly include the effect of water vapour, CO2, O3 and clouds:


The integral transmission functions for water vapour, CO2, O3 and clouds are given by

(A110)Dwv(z1,z2)=1+a1wv(β0Mwv(z1,z2))β1wv+a2wv(β0Mwv(z1,z2))β2wv+a3wv(β0Mwv(z1,z2))β3wv-1,(A111)DCO2(z1,z2)=1-0.1MCO2(z1,z2)10002×1+a0CO2a1CO2(β0MCO2(z1,z2))β1CO21+a0CO2(β0MCO2(z1,z2))β2CO2,(A112)DO3(z1,z2)=1-aO3(MO3(z1,z2))βO3,(A113)Dcld(z1,z2)=exp-|z-z|Hcldtop-Hcldbaseτcldonly inside cloud layers.

The first term on the right-hand side of Eq. (A111) ensures that the radiative forcing of CO2 increases with increasing CO2, in accordance with Hansen (2005) and Colman and McAvaney (2009). The effective absorber mass of water vapour, CO2 and O3 is computed from


The effect of the well-mixed greenhouse gases CH4, N2O and CFCs on longwave radiation is accounted for through a CO2 equivalent following Etminan et al. (2016):

(A117) CO 2 eq = CO 2 0 exp R ( CO 2 ) + R ( CH 4 ) + R ( N 2 O ) + R ( CFC 11 ) + R ( CFC 12 ) a 1 CO 2 - CO 2 0 2 + b 1 CO 2 - CO 2 0 + c 1 N 2 O + 5.36

The radiative forcing R of CO2, CH4 and N2O is computed as in Table 1 of Etminan et al. (2016) and the radiative forcing of CFC11 and CFC12 as in Table 3 of Myhre et al. (1998).

Similarly to shortwave radiation, longwave radiation is computed for each macro surface type individually. Atmospheric characteristics for each type within a given grid cell are the same but surface elevation and surface temperature are different. By default, for longwave radiation the atmosphere is divided into 15 vertical levels: 6 between surface and base of clouds, 3 between the base and top of clouds, 3 between the top of clouds and the tropopause and 3 in the stratosphere. Total longwave radiation fluxes in each grid cell are taken as cloud fraction weighted sum of clear sky and cloudy sky fluxes.

Longwave radiation parameters (in particular those related to the integral transmission functions) are tuned in offline mode using prescribed 3-D temperature and water vapour fields and cloud cover fraction from ERA-Interim reanalysis (Dee et al.2011), cloud top height and cloud optical thickness from the ISCCP climatology Rossow and Schiffer (1999), and O3 concentration from historical simulations of selected CMIP6 models. The target is to minimize the root mean square error in top-of-the-atmosphere and surface longwave radiation fluxes for clear sky and total sky. The resulting zonal mean net top-of-the-atmosphere longwave radiation is shown in Figs. A4 and A5, for clear-sky and total sky, respectively. Additional constraints on model parameters are provided by the longwave radiative kernels at the TOA and at the surface, which are available for several climate models (Shell et al.2008; Block and Mauritsen2013; Smith et al.2018; Pendergrass et al.2018). The TOA kernels computed using the CLIMBER-X longwave radiation code are compared to those available from the different climate models in Fig. A6, demonstrating that the simplified longwave radiation scheme employed in CLIMBER-X also represents well the sensitivities of longwave radiation fluxes to changes in temperature, humidity and CO2. The total greenhouse effect separation presented in Schmidt et al. (2010) is used as a further constraint on longwave radiation parameters. The effect of removing different longwave absorbers on the net longwave radiation at TOA is computed and compared to the values of Schmidt et al. (2010) in Table A8.

Table A7Longwave radiation parameters.

Download Print Version | Download XLSX

Schmidt et al. (2010)

Table A8Effect of removing different longwave absorbers on the net longwave absorbed by the circa 1980 atmosphere. The unit is W m−2. The reference values are from Schmidt et al. (2010).

Download Print Version | Download XLSX

Figure A4Zonal mean net longwave radiation at the top of the atmosphere for clear-sky conditions as computed by the CLIMBER-X radiation code in an offline set-up compared to ERA-Interim reanalysis (Dee et al.2011) and CERES satellite observations (Loeb et al.2018). Panels (a) and (b) are for December–January–February and panels (c) and (d) are for June–July–August. Panels (a) and (c) show absolute values, while panels (b) and (d) show model biases.

Figure A5Zonal mean net longwave radiation at the top of the atmosphere as computed by the CLIMBER-X radiation code in an offline set-up compared to ERA-Interim reanalysis (Dee et al.2011) and CERES satellite observations (Loeb et al.2018). Panels (a) and (b) are for December–January–February and panels (c) and (d) are for June–July–August. Panels (a) and (c) show absolute values, while panels (b) and (d) show model biases.

Figure A6Longwave radiative kernels at the top of the atmosphere for (a) surface temperature warming by 1 K, (b) atmospheric CO2 doubling, (c) air temperature warming by 1 K and (d) water vapour increase corresponding to 1 K warming in the atmosphere. The kernels computed using the CLIMBER-X radiation code are compared to kernels available from different climate models (Shell et al.2008; Block and Mauritsen2013; Pendergrass et al.2018). Solid lines are for total-sky conditions, while dashed lines are for clear-sky conditions.

Appendix B: GOLDSTEIN detailed description

The horizontal ocean velocity in GOLDSTEIN is diagnosed from a frictional–geostrophic balance:


The drag coefficient μ is composed of a uniform background value μ0, which, for barotropic velocities, is enhanced by a factor of 3 in waters shallower than 1000 m. The definition of the wind or sea ice stress τ is given in Appendix A5 above. Hydrostatic balance is assumed:

(B3) p z = - g ρ ( θ , S , p 0 ( z ) ) ,

where p0(z)=0.1z. Seawater density, ρ(θ,S,p0(z)), is computed using the UNESCO equation of state of Millero and Poisson (1981) with the bulk secant modulus expressed in terms of potential temperature using the coefficients of Jackett and McDougall (1995). For computational efficiency only 7 of the originally 26 terms are kept in the bulk secant modulus polynomial. A comparison between the full equation of state of Jackett and McDougall (1995) and the truncated version showed negligible differences in the simulated ocean state. The vertical velocity is then derived from the continuity equation:

(B4) u = 0 .

Equations (B1) through (B4) are solved by separation of the velocity into a barotropic and a baroclinic component, as described in detail in Edwards et al. (1998) and Müller et al. (2006). More details on the treatment of barotropic flow around islands are given in Edwards and Shepherd (2002).

The transport equation for tracers takes the following form:

(B5) X t + ( u X ) = ( A X ) + SMS .

X is the tracer concentration, SMS represents source minus sinks and A is the diffusive mixing tensor:

(B6) A = K I 0 ( K I - κ ) S λ 0 K I ( K I - κ ) S ϕ ( K I + κ ) S λ ( K I + κ ) S ϕ K I S 2 + K D ,

with KI and KD the isopycnal and diapycnal diffusivities, respectively, κ the Gent–McWilliams parameter which parameterizes sub-scale eddies and the slope of the isopycnals is:

(B7) S = S λ , S ϕ , 0 = - 1 R e cos ϕ ρ λ ρ z , - 1 R e ρ ϕ ρ z , 0 .

Diapycnal diffusivity is a prescribed function of depth:

(B8) K D = K D min + arctan z - z ref 1000 - arctan - z ref 1000 arctan 5000 - z ref 1000 - arctan - z ref 1000 K D max - K D min .

If the stratification of a water column is statically unstable, convective adjustment is applied using the scheme of Rahmstorf (1993). A simple mixed-layer scheme based on Kraus and Turner (1967) is used in the model.

The equations are discretized on a staggered Arakawa C grid.

Table B1Ocean model parameters.

Download Print Version | Download XLSX

Appendix C: SISIM detailed description

Sea ice thermodynamics is based on the zero-layer model of Semtner (1976). The surface energy fluxes are computed separately over sea ice and over open ocean. Over sea ice, the surface energy balance equation is written as

(C1) ( 1 - α ) SW + ϵ LW - LW - SH - LE - G = 0 ,

where α is surface albedo, SW is the incoming shortwave radiation, ϵ is the surface emissivity for longwave radiation, LW and LW are the incoming and outgoing longwave radiation at the surface, SH is the sensible heat flux, LE is the latent heat flux and G the heat flux into the snow/ice. Equation (C1) is then solved for the skin temperature, T, using the formulations for the energy fluxes described next. Bare sea ice albedo is temperature-dependent following


The snow albedo scheme is the same as used in the land model and includes a dependence on snow grain size and dust and soot concentration following Dang et al. (2015). The fraction of sea ice covered by snow is computed from snow thickness as

(C4) f snow = h snow h snow + 0.02 .

The surface emitted longwave radiation is given by the Stefan–Boltzmann law:

(C5) LW = ϵ σ T 4 .

The sensible heat flux is computed from the temperature gradient between the skin and near-surface air, using the bulk aerodynamic formula:

(C6) SH = ρ a c p r aer ( T - T 2 m ) ,

where ρa is air density, cp is the specific heat of air, raer is the aerodynamic resistance and T2 m is the near-surface air temperature. Similarly, the latent heat flux over sea ice is expressed in terms of the specific humidity gradient between the surface and near-surface air:

(C7) LE = L ρ a r aer q sat ( T ) - q 2 m .

L is the latent heat of sublimation, qsat is the specific humidity at saturation and q2 m is the specific humidity of near-surface air. The aerodynamic resistance is computed from wind speed, surface exchange coefficient and bulk Richardson number following Willeit and Ganopolski (2016).

The conductive heat flux into the snow/ice (G) is computed as

(C8) G = λ eff T - T f ,

where Tf is the salinity-dependent freezing temperature at the base of the ice (Millero and Poisson1981):

(C9) T f = - 0.0575 S o + 0.0017 S o 1.5 - 0.0002 S o 2 + T 0

and the effective snow/ice heat conductivity is

(C10) λ eff = k λ λ snow h snow + h ice λ snow / λ ice .

The kλ factor accounts for sub-grid ice thickness distribution following Eq. (2) in (Fichefet and Maqueda1997). The prognostic terms in T in the formulation of the surface energy fluxes are then linearized using Taylor series expansion assuming that the temperature at the new time step, T,n+1=T,n+ΔT with ΔTT:


Equation (C1) can then be solved explicitly for the skin temperature at the new time step, T,n+1. If the skin temperature is above freezing the surface energy fluxes are diagnosed first with the skin temperature greater then 0 C and then with skin temperature set to 0 C. The difference between the sum of the energy fluxes is then used to melt snow and/or ice.

Whether bottom ice accretion or ablation occurs depends on the sign of the net energy flux at the ice–water interface, which is determined by the balance between the conductive heat flux through the snow/ice layer (G) and the turbulent heat flux between the ice and the seawater below (McPhee1992; Weaver et al.2001):

(C13) H = C h u ρ w c w T f - T o .

Ch is a constant exchange coefficient, u the (constant) friction velocity, ρw water density and cw the specific heat capacity of water.

Over seawater, the heat flux into the ocean is derived as the residual of the radiation and surface energy fluxes, which are computed similarly to the surface energy fluxes over sea ice, but using sea surface temperature instead of skin temperature. The albedo of seawater for diffuse radiation is 0.06, while the clear-sky albedo includes a dependence on the cosine of the solar zenith angle (μ):

(C14) α = min 0.2 , 0.03 μ .

The longwave emissivity of water is set to 0.98. Sea ice in leads forms whenever the top-layer ocean temperature drops below the freezing point of seawater.

Table C1Sea ice model parameters.

Download Print Version | Download XLSX

Sub-grid-scale thermodynamic processes of sea ice growth and melt are assumed to affect the sea ice concentration (fice) within a grid cell following Marsland et al. (2003) and Fichefet and Maqueda (1997):

(C15) f ice t = Δ f ice | frz + Δ f ice | melt .

When freezing occurs over open water the sea ice concentration increases at a rate given by

(C16) Δ f ice | frz = 1 - f ice 2 1 - f ice Δ h ice | frz h 0 ,

where Δhice|frz is the thickness of new sea ice formed and h0 is an arbitrary demarcation thickness. Melting of sea ice results in a decrease in sea ice concentration:

(C17) Δ f ice | melt = f ice Δ h ice | melt 2 h ice ,

where Δhice|melt is the change in sea ice thickness due to melting. This formulation is based on the assumption that sea ice thickness within a grid cell has a uniform distribution between 0 and 2hice.

Sea ice drift velocities are computed from the momentum balance equation (e.g. Hibler1979):

(C18) m u i t = σ + τ a + τ o - k ^ × m f u i - m g H o ,

where ui=ui,vi is the sea ice velocity vector, m is the mass of ice and snow per unit area, f is the Coriolis parameter, Ho is the sea surface height and τa and τo are the wind and ocean stresses, respectively. The mechanical properties of the ice are represented by the internal stress tensor σ, and the elastic–viscous–plastic rheology of Hunke and Dukowicz (1997); Bouillon et al. (2009) is used as a constitutive law relating the internal stress to the strain rate. The numerical solution follows the implementation in the GFDL model SIS2 (Adcroft et al.2019; Delworth et al.2006). The sea surface height is diagnosed to a first approximation from the seawater density above a reference depth of 1500 m. The zonal and meridional components of the wind stress on sea ice are computed as


with the symbols defined in Appendix A2 above. The components of the stress exerted by the ocean on sea ice are given by


where Cdw is the drag coefficient between the sea ice and water and (uo,vo) the top-layer ocean velocity. The derived sea ice velocities are then used to advect the grid cell mean ice and snow thicknesses and the sea ice concentration using a flux-corrected transport scheme (Zalesak1979). No explicit diffusion is applied.

Appendix D: PALADYN surface albedo

Several changes to the surface albedo scheme have been introduced into the model compared to the description in Willeit and Ganopolski (2016).

The parameterization of sub-grid snow cover fraction has been updated considering also the effect of topographic roughness following Niu and Yang (2007) and Roesch et al. (2001):

(D1) f sn = tanh h sn 10 z 0 h sn h sn + 1 × 10 - 4 σ oro ,

where hsn is snow thickness, z0 is the surface roughness length and σoro is the sub-grid standard deviation of the orography.

Figure D1Snow grain size parameterization as a function of skin temperature and snowfall rate.


The snow albedo scheme has been refined with the inclusion of the effect of dust and soot on snow albedo following Dang et al. (2015) and the snow grain size parameterization has been retuned using output of the regional climate model MARv3.6 simulations for Greenland (Fettweis et al.2017). The resulting dependence of snow grain size on climatological values of skin temperature and snowfall rate are shown in Fig. D1.

Code and data availability

The source code of CLIMBER-X v1.0 is archived on Zenodo (, Willeit2022), together with the input data and the instructions to install and run the model. The FFTW library, which is already included in CLIMBER-X, is available from (Frigo and Johnson2005). CMIP6 model data are licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (, last access: 8 December 2021) and can be accessed through the ESGF nodes (for instance, last access: 10 November 2021). ERA-Interim data are available from (Dee et al.2011), while ERA5 reanalysis data can be downloaded from the Copernicus Climate Data Store at (Hersbach et al.2019).

Author contributions

AG and MW conceived the model. AG developed the atmosphere model SESAM, with contributions from MW. NRE provided the GOLDSTEIN model code and supported the implementation of the ocean model. MW developed the sea ice model SISIM. MW coupled the different model components and tuned and tested the model. AR contributed to the model design and provided tools for model I/O, mapping/interpolation and ensemble generation. MW performed the model simulations, prepared the figures and wrote the paper, with contributions from all the authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Matteo Willeit was supported by the German Science Foundation (DFG) grant GA 1202/2-1 and by the BMBF-funded project PalMod. Alexander Robinson was funded by the Ramón y Cajal Programme of the Spanish Ministry for Science, Innovation and Universities (grant no. RYC-2016-20587). We thank the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP5 and CMIP6. We thank the climate modelling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP5, CMIP6 and ESGF. Xavier Fettweiss is thanked for providing the grain size data from simulations of the regional climate model MAR. Data from the RAPID AMOC monitoring project are funded by the Natural Environment Research Council and are freely available from (last access: 15 January 2022​​​​​​​).

Financial support

This research has been supported by the Bundesministerium für Bildung und Forschung (grant no. PalMod), the Deutsche Forschungsgemeinschaft (grant no. GA 1202/2-1), and the Ministerio de Ciencia e Innovación (grant no. RYC-2016-20587).

The publication of this article was funded by the Open Access Fund of the Leibniz Association.

Review statement

This paper was edited by Riccardo Farneti and reviewed by three anonymous referees.


Adcroft, A., Anderson, W., Balaji, V., Blanton, C., Bushuk, M., Dufour, C. O., Dunne, J. P., Griffies, S. M., Hallberg, R., Harrison, M. J., Held, I. M., Jansen, M. F., John, J. G., Krasting, J. P., Langenhorst, A. R., Legg, S., Liang, Z., McHugh, C., Radhakrishnan, A., Reichl, B. G., Rosati, T., Samuels, B. L., Shao, A., Stouffer, R., Winton, M., Wittenberg, A. T., Xiang, B., Zadeh, N., and Zhang, R.: The GFDL Global Ocean and Sea Ice Model OM4.0: Model Description and Simulation Features, J. Adv. Model. Earth Sy., 11, 3167–3211,, 2019. a, b

Adkins, J. F., McIntyre, K., and Schrag, D. P.: The salinity, temperature, and δ18O of the glacial deep ocean, Science, 298, 1769–1773,, 2002. a

Bala, G., Caldeira, K., Mirin, A., Wickett, M., Delire, C., and Phillips, T. J.: Biogeophysical effects of CO2 fertilization on global climate, Tellus B, 58, 620–627,, 2006. a

Bauer, E. and Ganopolski, A.: Aeolian dust modeling over the past four glacial cycles with CLIMBER-2, Global Planet. Change, 74, 49–60,, 2010. a, b

Bauer, E., Petoukhov, V., Ganopolski, A., and Eliseev, A. V.: Climatic response to anthropogenic sulphate aerosols versus well-mixed greenhouse gases from 1850 to 2000 AD in CLIMBER-2, Tellus B, 60B, 82–97,, 2008. a, b

Bereiter, B., Shackleton, S., Baggenstos, D., Kawamura, K., and Severinghaus, J.: Mean global ocean temperatures during the last glacial transition, Nature, 553, 39–44,, 2018. a

Block, K. and Mauritsen, T.: Forcing and feedback in the MPI-ESM-LR coupled model under abruptly quadrupled CO2, J. Adv. Model. Earth Sy., 5, 676–691,, 2013. a, b, c, d, e

Bock, L., Lauer, A., Schlund, M., Barreiro, M., Bellouin, N., Jones, C., Meehl, G. A., Predoi, V., Roberts, M. J., and Eyring, V.: Quantifying Progress Across Different CMIP Phases With the ESMValTool, J. Geophys. Res.-Atmos., 125, e2019JD032321​​​​​​​,, 2020. a

Bohm, E., Lippold, J., Gutjahr, M., Frank, M., Blaser, P., Antz, B., Fohlmeister, J., Frank, N., Andersen, M. B., and Deininger, M.: Strong and deep Atlantic meridional overturning circulation during the last glacial cycle, Nature, 517, 73–76,, 2015. a

Bonan, G. B.: Forests and climate change: forcings, feedbacks, and the climate benefits of forests, Science, 320, 1444–1449,, 2008. a

Bony, S., Colman, R., Kattsov, V. M., Allan, R. P., Bretherton, C. S., Dufresne, J.-L. L., Hall, A., Hallegatte, S., Holland, M. M., Ingram, W., Randall, D. a., Soden, B. J., Tselioudis, G., and Webb, M. J.: How Well Do We Understand and Evaluate Climate Change Feedback Processes?, J. Climate, 19, 3445–3482,, 2006. a, b

Bouillon, S., Morales Maqueda, M. Á., Legat, V., and Fichefet, T.: An elastic–viscous–plastic sea ice model formulated on Arakawa B and C grids, Ocean Model., 27, 174–184,, 2009. a, b

Brown, J., Ferrians, O., Heginbottom, J. A., and Melnikov, E.: Circum-Arctic Map of Permafrost and Ground-Ice Conditions, Version 2, Boulder, Colorado USA, NSIDC: National Snow and Ice Data Center [data set],, 1998. a

Bryan, K. and Lewis, L. J.: A water mass model of the World Ocean, J. Geophys. Res., 84, 2503–2517​​​​​​​,, 1979. a

Burton, C., Betts, R., Cardoso, M., Feldpausch, T. R., Harper, A., Jones, C. D., Kelley, D. I., Robertson, E., and Wiltshire, A.: Representation of fire, land-use change and vegetation dynamics in the Joint UK Land Environment Simulator vn4.9 (JULES), Geosci. Model Dev., 12, 179–193,, 2019. a

Caballero, R. and Hanley, J.: Midlatitude eddies, storm-track diffusivity, and poleward moisture transport in warm climates, J. Atmos. Sci., 69, 3237–3250,, 2012. a

Calov, R., Ganopolski, A., Claussen, M., Petoukhov, V., and Greve, R.: Transient simulation of the last glacial inception. Part I: glacial inception as a bifurcation in the climate system, Clim. Dynam., 24, 545–561,, 2005. a

Charney, J., Arakawa, A., Baker, D., Bolin, B., Dickinson, R., Goody, R., Leith, C., Stommel, H., and Wunsch, C.: Carbon Dioxide and Climate: A Scientific Assessment, Tech. Rep., National Academy of Sciences, Washington, D.C.,, 1979. a

Charney, J. G. and Eliassen, A.: A Numerical Method for Predicting the Perturbations of the Middle Latitude Westerlies, Tellus, 1, 38–54,, 1949. a, b

Claussen, M., Mysak, L., Weaver, A., Crucifix, M., Fichefet, T., Loutre, M. F., Weber, S., Alcamo, J., Alexeev, V., Berger, A., Calov, R., Ganopolski, A., Goosse, H., Lohmann, G., Lunkeit, F., Mokhov, I., Petoukhov, V., Stone, P., and Wang, Z.: Earth system models of intermediate complexity: Closing the gap in the spectrum of climate system models, Clim. Dynam., 18, 579–586,, 2002. a

Colman, R. and McAvaney, B.: Climate feedbacks under a very broad range of forcing, Geophys. Res. Lett., 36, 1–5​​​​​​​,, 2009. a, b

Colman, R., Fraser, J., and Rotstayn, L.: Climate feedbacks in a general circulation model incorporating prognostic clouds, Clim. Dynam., 18, 103–122,, 2001. a

Crook, J. A., Forster, P. M., and Stuber, N.: Spatial patterns of modeled climate feedback and contributions to temperature response and polar amplification, J. Climate, 24, 3575–3592,, 2011. a

Dang, C., Brandt, R. E., and Warren, S. G.: Parameterizations for narrowband and broadband albedo of pure snow and snow containing mineral dust and black carbon, J. Geophys. Res.-Atmos., 120, 5446–5468,, 2015. a, b, c, d

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597,, 2011. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Delworth, T. L., Broccoli, A. J., Rosati, A., Stouffer, R. J., Balaji, V., Beesley, J. A., Cooke, W. F., Dixon, K. W., Dunne, J., Dunne, K. A., Durachta, J. W., Findell, K. L., Ginoux, P., Gnanadesikan, A., Gordon, C. T., Griffies, S. M., Gudgel, R., Harrison, M. J., Held, I. M., Hemler, R. S., Horowitz, L. W., Klein, S. A., Knutson, T. R., Kushner, P. J., Langenhorst, A. R., Lee, H.-C., Lin, S.-J., Lu, J., Malyshev, S. L., Milly, P. C. D., Ramaswamy, V., Russell, J., Schwarzkopf, M. D., Shevliakova, E., Sirutis, J. J., Spelman, M. J., Stern, W. F., Winton, M., Wittenberg, A. T., Wyman, B., Zeng, F., and Zhang, R.: GFDL's CM2 Global Coupled Climate Models. Part I: Formulation and Simulation Characteristics, J. Climate, 19, 643–674,, 2006. a, b

Durack, P. J., Wijffels, S. E., and Matear, R. J.: Ocean Salinities Reveal Strong Global Water Cycle Intensification During 1950 to 2000, Science, 336, 455–458,, 2012. a, b

Eby, M., Weaver, A. J., Alexander, K., Zickfeld, K., Abe-Ouchi, A., Cimatoribus, A. A., Crespin, E., Drijfhout, S. S., Edwards, N. R., Eliseev, A. V., Feulner, G., Fichefet, T., Forest, C. E., Goosse, H., Holden, P. B., Joos, F., Kawamiya, M., Kicklighter, D., Kienert, H., Matsumoto, K., Mokhov, I. I., Monier, E., Olsen, S. M., Pedersen, J. O. P., Perrette, M., Philippon-Berthier, G., Ridgwell, A., Schlosser, A., Schneider von Deimling, T., Shaffer, G., Smith, R. S., Spahni, R., Sokolov, A. P., Steinacher, M., Tachiiri, K., Tokos, K., Yoshimori, M., Zeng, N., and Zhao, F.: Historical and idealized climate model experiments: an intercomparison of Earth system models of intermediate complexity, Clim. Past, 9, 1111–1140,, 2013. a

ECCO Consortium, Fukumori, I., Wang, O., Fenty, I., Forget, G., Heimbach, P., and Ponte, R. M.: ECCO Ocean Mixed Layer Depth – Monthly Mean 0.5 Degree (Version 4 Release 4), Ver. V4r4, NASA [data set],, 2021. a

Edwards, N. and Shepherd, J.: Bifurcations of the thermohaline circulation in a simplified three-dimensional model of the world ocean and the effects of inter-basin connectivity, Clim. Dynam., 19, 31–42,, 2002. a, b

Edwards, N. R. and Marsh, R.: Uncertainties due to transport-parameter sensitivity in an efficient 3-D ocean-climate model, Clim. Dynam., 24, 415–433,, 2005. a, b

Edwards, N. R., Willmott, A. J., and Killworth, P. D.: On the Role of Topography and Wind Stress on the Stability of the Thermohaline Circulation, J. Phys. Oceanogr., 28, 756–778,<0756:OTROTA>2.0.CO;2, 1998. a, b

Etminan, M., Myhre, G., Highwood, E. J., and Shine, K. P.: Radiative forcing of carbon dioxide, methane, and nitrous oxide: A significant revision of the methane radiative forcing, Geophys. Res. Lett., 43, 12614–12623,, 2016. a, b, c

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,, 2016. a

Falloon, P. D., Dankers, R., Betts, R. A., Jones, C. D., Booth, B. B. B., and Lambert, F. H.: Role of vegetation change in future climate under the A1B scenario and a climate stabilisation scenario, using the HadCM3C Earth system model, Biogeosciences, 9, 4739–4756,, 2012. a

Farneti, R. and Vallis, G. K.: An Intermediate Complexity Climate Model (ICCMp1) based on the GFDL flexible modelling system, Geosci. Model Dev., 2, 73–88,, 2009. a

Fasullo, J. T. and Trenberth, K. E.: The annual cycle of the energy budget. Part II: Meridional structures and poleward transports, J. Climate, 21, 2313–2325,, 2008. a, b, c

Fedorovich, E. and Shapiro, A.: Structure of numerically simulated katabatic and anabatic flows along steep slopes, Acta Geophys., 57, 981–1010,, 2009. a

Feigelson, E., Ginzburg, A., Krasnokutskaya, L., and Petoukhov, V.: Effects of clouds on the radiative heat exchange in the atmosphere, Geofís. Int., 15, 293–326,, 1975. a, b

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033,, 2017. a, b

Fichefet, T. and Maqueda, M. A. M.: Sensitivity of a global sea ice model to the treatment of ice thermodynamics and dynamics, J. Geophys. Res.-Oceans, 102, 12609–12646,, 1997. a, b, c

Fraedrich, K., Kirk, E., Luksch, U., and Lunkeit, F.: The portable university model of the atmosphere (PUMA): Storm track dynamics and low-frequency variability, Meteorol. Z., 14, 735–745,, 2005. a

Frajka-Williams, E., Moat, B., Smeed, D., Rayner, D., Johns, W., Baringer, M., Volkov, D., and Collins, J.: Atlantic meridional overturning circulation observed by the RAPID-MOCHA-WBTS (RAPID-Meridional Overturning Circulation and Heatflux Array-Western Boundary Time Series) array at 26N from 2004 to 2020 (v2020.1), National Oceanography Centre [data set],, 2021. a, b

Frierson, D. M., Lu, J., and Chen, G.: Width of the Hadley cell in simple and comprehensive general circulation models, Geophys. Res. Lett., 34, 1–5​​​​​​​,, 2007. a

Frigo, M. and Johnson, S.: The Design and Implementation of FFTW3, P. IEEE, 93, 216–231, 2005. a, b

Ganopolski, A. and Brovkin, V.: Simulation of climate, ice sheets and CO2 evolution during the last four glacial cycles with an Earth system model of intermediate complexity, Clim. Past, 13, 1695–1716,, 2017. a

Ganopolski, A., Rahmstorf, S., Petoukhov, V., and Claussen, M.: Simulation of modern and glacial climates with a coupled global model of intermediate complexity, Nature, 391, 351–356,, 1998. a

Ganopolski, A., Petoukhov, V., Rahmstorf, S., Brovkin, V., Claussen, M., Eliseev, A., and Kubatzki, C.: CLIMBER-2: a climate system model of intermediate complexity. Part II: model sensitivity, Clim. Dynam., 17, 735–751,, 2001. a

Ganopolski, A., Winkelmann, R., and Schellnhuber, H. J.: Critical insolation–CO2 relation for diagnosing past and future glacial inception, Nature, 529, 200–203,, 2016. a

Gent, P. R. and Mcwilliams, J. C.: Isopycnal Mixing in Ocean Circulation Models, J. Phys. Oceanogr., 20, 150–155,<0150:IMIOCM>2.0.CO;2, 1990. a

Gerdes, R., Köberle, C., and Willebrand, J.: The influence of numerical advection schemes on the results of ocean general circulation models, Clim. Dynam., 5, 211–226,, 1991. a

Goosse, H., Brovkin, V., Fichefet, T., Haarsma, R., Huybrechts, P., Jongma, J., Mouchet, A., Selten, F., Barriat, P.-Y., Campin, J.-M., Deleersnijder, E., Driesschaert, E., Goelzer, H., Janssens, I., Loutre, M.-F., Morales Maqueda, M. A., Opsteegh, T., Mathieu, P.-P., Munhoven, G., Pettersson, E. J., Renssen, H., Roche, D. M., Schaeffer, M., Tartinville, B., Timmermann, A., and Weber, S. L.: Description of the Earth system model of intermediate complexity LOVECLIM version 1.2, Geosci. Model Dev., 3, 603–633,, 2010. a

Greve, R.: Application of a Polythermal Three-Dimensional Ice Sheet Model to the Greenland Ice Sheet: Response to Steady-State and Transient Climate Scenarios, J. Climate, 10, 901–918,<0901:AOAPTD>2.0.CO;2, 1997. a

Griffies, S. M.: The Gent–McWilliams Skew Flux, J. Phys. Oceanogr., 28, 831–841,<0831:TGMSF>2.0.CO;2, 1998. a

Hansen, J.: Efficacy of climate forcings, J. Geophys. Res., 110, D18104,, 2005. a

Hansen, J., Russell, G., Rind, D., Stone, P., Lacis, A., Lebedeff, S., Ruedy, R., and Travis, L.: Efficient Three-Dimensional Global Models for Climate Studies: Models I and II, Mon. Weather Rev., 111, 609–662,<0609:ETDGMF>2.0.CO;2, 1983. a, b

Hawkins, E., Smith, R. S., Allison, L. C., Gregory, J. M., Woollings, T. J., Pohlmann, H., and De Cuevas, B.: Bistability of the Atlantic overturning circulation in a global climate model and links to ocean freshwater transport, Geophys. Res. Lett., 38, 1–6​​​​​​​,, 2011. a

Held, I. M.: Stationary and quasi-stationary eddies in the extratropical troposphere: Theory, in: Large-Scale Dynamical Processes in the Atmosphere, edited by: Hoskins, B. and Pearce, R. P., Academic Press, 127–168, ISBN-10 0123566800, ISBN-13 978-0123566805, 1983. a

Held, I. M. and Soden, B. J.: Robust Responses of the Hydrological Cycle to Global Warming, J. Climate, 19, 5686–5699,, 2006. a

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., and Thépaut, J.-N.: ERA5 monthly averaged data on single levels from 1959 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set],, 2019. a

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J. N.: The ERA5 global reanalysis, Q. J. Roy. Meteorol. Soc., 146, 1999–2049,, 2020. a

Hibler, W. D.: A Dynamic Thermodynamic Sea Ice Model, J. Phys. Oceanogr., 9, 815–846,<0815:ADTSIM>2.0.CO;2, 1979. a, b

Holden, P. B., Edwards, N. R., Fraedrich, K., Kirk, E., Lunkeit, F., and Zhu, X.: PLASIM–GENIE v1.0: a new intermediate complexity AOGCM, Geosci. Model Dev., 9, 3347–3361,, 2016. a, b

Holton, J. R.: Chapter 7 Atmospheric oscillations: Linear perturbation theory, in: An Introduction to Dynamic Meteorology, edited by: Holton, J. R.​​​​​​​, vol. 88, Academic Press, 182–227,, 2004. a

Hoskins, B. J. and Valdes, P. J.: On the Existence of Storm-Tracks, J. Atmos. Sci., 47, 1854–1864,<1854:OTEOST>2.0.CO;2, 1990. a, b

Hu, Y., Huang, H., and Zhou, C.: Widening and weakening of the Hadley circulation under global warming, Sci. Bull., 63, 640–644,, 2018. a

Hunke, E. C. and Dukowicz, J. K.: An elastic-viscous-plastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867,<1849:AEVPMF>2.0.CO;2, 1997. a, b

Ilyina, T., Six, K. D., Segschneider, J., Maier-Reimer, E., Li, H., and Núñez-Riboni, I.: Global ocean biogeochemistry model HAMOCC: Model architecture and performance as component of the MPI-Earth system model in different CMIP5 experimental realizations, J. Adv. Model. Earth Sy., 5, 287–315,, 2013. a

IPCC: Annex II: Climate System Scenario Tables, edited by: Prather, M., Flato, G., Friedlingstein, P., Jones, C., Lamarque, J.-F., Liao, H., and Rasch, P., in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1395–1446,, 2013. a

IPCC: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, in press,, 2021. a

Jackett, D. R. and McDougall, T. J.: Minimal Adjustment of Hydrographic Profiles to Achieve Static Stability, J. Atmos. Ocean. Tech., 12, 381–389,<0381:maohpt>;2, 1995. a, b

Jackson, L. C., Kahana, R., Graham, T., Ringer, M. A., Woollings, T., Mecking, J. V., and Wood, R. A.: Global and European climate impacts of a slowdown of the AMOC in a high resolution GCM, Clim. Dynam., 45, 3299–3316,, 2015. a

Johns, W. E., Baringer, M. O., Beal, L. M., Cunningham, S. A., Kanzow, T., Bryden, H. L., Hirschi, J. J., Marotzke, J., Meinen, C. S., Shaw, B., and Curry, R.: Continuous, array-based estimates of atlantic ocean heat transport at 26.5 N, J. Climate, 24, 2429–2449,, 2011. a

Kageyama, M., Albani, S., Braconnot, P., Harrison, S. P., Hopcroft, P. O., Ivanovic, R. F., Lambert, F., Marti, O., Peltier, W. R., Peterschmitt, J.-Y., Roche, D. M., Tarasov, L., Zhang, X., Brady, E. C., Haywood, A. M., LeGrande, A. N., Lunt, D. J., Mahowald, N. M., Mikolajewicz, U., Nisancioglu, K. H., Otto-Bliesner, B. L., Renssen, H., Tomas, R. A., Zhang, Q., Abe-Ouchi, A., Bartlein, P. J., Cao, J., Li, Q., Lohmann, G., Ohgaito, R., Shi, X., Volodin, E., Yoshida, K., Zhang, X., and Zheng, W.: The PMIP4 contribution to CMIP6 – Part 4: Scientific objectives and experimental design of the PMIP4-CMIP6 Last Glacial Maximum experiments and PMIP4 sensitivity experiments, Geosci. Model Dev., 10, 4035–4055,, 2017. a

Kageyama, M., Harrison, S. P., Kapsch, M.-L., Lofverstrom, M., Lora, J. M., Mikolajewicz, U., Sherriff-Tadano, S., Vadsaria, T., Abe-Ouchi, A., Bouttes, N., Chandan, D., Gregoire, L. J., Ivanovic, R. F., Izumi, K., LeGrande, A. N., Lhardy, F., Lohmann, G., Morozova, P. A., Ohgaito, R., Paul, A., Peltier, W. R., Poulsen, C. J., Quiquet, A., Roche, D. M., Shi, X., Tierney, J. E., Valdes, P. J., Volodin, E., and Zhu, J.: The PMIP4 Last Glacial Maximum experiments: preliminary results and comparison with the PMIP3 simulations, Clim. Past, 17, 1065–1089,, 2021. a, b, c

Klemann, V., Martinec, Z., and Ivins, E. R.: Glacial isostasy and plate motion, J. Geodyn., 46, 95–103,, 2008. a

Köhler, P., Nehrbass-Ahles, C., Schmitt, J., Stocker, T. F., and Fischer, H.: A 156 kyr smoothed history of the atmospheric greenhouse gases CO2, CH4, and N2O and their radiative forcing, Earth Syst. Sci. Data, 9, 363–387,, 2017. a

Krapp, M., Robinson, A., and Ganopolski, A.: SEMIC: an efficient surface energy and mass balance model applied to the Greenland ice sheet, The Cryosphere, 11, 1519–1535,, 2017. a, b

Kraus, E. B. and Turner, J. S.: A one-dimensional model of the seasonal thermocline II. The general theory and its consequences, Tellus, 19, 98–106,, 1967. a, b

Lacis, A. A. and Hansen, J.: A Parameterization for the Absorption of Solar Radiation in the Earth's Atmosphere, J. Atmos. Sci., 31, 118–133,<0118:APFTAO>2.0.CO;2, 1974. a

Lenton, T. M., Marsh, R., Price, A. R., Lunt, D. J., Aksenov, Y., Annan, J. D., Cooper-Chadwick, T., Cox, S. J., Edwards, N. R., Goswami, S., Hargreaves, J. C., Harris, P. P., Jiao, Z., Livina, V. N., Payne, A. J., Rutt, I. C., Shepherd, J. G., Valdes, P. J., Williams, G., Williamson, M. S., and Yool, A.: Effects of atmospheric dynamics and ocean resolution on bi-stability of the thermohaline circulation examined using the Grid ENabled Integrated Earth system modelling (GENIE) framework, Clim. Dynam., 29, 591–613,, 2007. a

Levis, S., Foley, J. A., and Pollard, D.: Potential high-latitude vegetation feedbacks on CO2-induced climate change, Geophys. Res. Lett., 26, 747–750,, 1999. a

Levitus, S., Antonov, J. I., Boyer, T. P., Baranova, O. K., Garcia, H. E., Locarnini, R. A., Mishonov, A. V., Reagan, J. R., Seidov, D., Yarosh, E. S., and Zweng, M. M.: World ocean heat content and thermosteric sea level change (0–2000 m), 1955–2010, Geophys. Res. Lett., 39, 1–5​​​​​​​,, 2012. a, b

Levitus, S., Boyer, T. P., and Garcia, Hernan E. Locarnini, Ricardo A. Zweng, Melissa M. Mishonov, Alexey V. Reagan, James R. Antonov, John I. Baranova, Olga K. Biddle, Mathew Hamilton, Melanie Johnson, Daphne R. Paver, Christopher R. Seidov, D.: World Ocean Atlas 2013 (NCEI Accession 0114815), NCEI [data set],, 2015. a, b

Lhardy, F., Bouttes, N., Roche, D. M., Crosta, X., Waelbroeck, C., and Paillard, D.: Impact of Southern Ocean surface conditions on deep ocean circulation during the LGM: a model analysis, Clim. Past, 17, 1139–1159,, 2021. a, b

Loeb, N. G., Doelling, D. R., Wang, H., Su, W., Nguyen, C., Corbett, J. G., Liang, L., Mitrescu, C., Rose, F. G., and Kato, S.: Clouds and the Earth'S Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) top-of-atmosphere (TOA) edition-4.0 data product, J. Climate, 31, 895–918,, 2018. a, b, c, d, e

Lucazeau, F.: Analysis and Mapping of an Updated Terrestrial Heat Flow Data Set, Geochem. Geophy. Geosy., 20, 4001–4024,, 2019. a

Ma, L., Hurtt, G. C., Chini, L. P., Sahajpal, R., Pongratz, J., Frolking, S., Stehfest, E., Klein Goldewijk, K., O'Leary, D., and Doelman, J. C.: Global rules for translating land-use change (LUH2) to land-cover change for CMIP6 using GLM2, Geosci. Model Dev., 13, 3203–3220,, 2020. a

Maier-Reimer, E. and Hasselmann, K.: Transport and storage of CO2 in the ocean – an inorganic ocean-circulation carbon cycle model, Clim. Dynam., 2, 63–90,, 1987. a

Manabe, S. and Stouffer, R. J.: Two Stable Equilibria of a Coupled Ocean-Atmosphere Model, J. Climate, 841–866,<0841:TSEOAC>2.0.CO;2​​​​​​​, 1988. a

Marsh, R., Müller, S. A., Yool, A., and Edwards, N. R.: Incorporation of the C-GOLDSTEIN efficient climate model into the GENIE framework: “eb_go_gs” configurations of GENIE, Geosci. Model Dev., 4, 957–992,, 2011. a, b

Marsland, S., Haak, H., Jungclaus, J., Latif, M., and Röske, F.: The Max-Planck-Institute global ocean/sea ice model with orthogonal curvilinear coordinates, Ocean Model., 5, 91–127,, 2003. a

Martinec, Z., Klemann, V., van der Wal, W., Riva, R. E., Spada, G., Sun, Y., Melini, D., Kachuck, S. B., Barletta, V., Simon, K., A, G., and James, T. S.: A benchmark study of numerical implementations of the sea level equation in GIA modelling, Geophys. J. Int., 215, 389–414,, 2018. a

Matthes, K., Funke, B., Andersson, M. E., Barnard, L., Beer, J., Charbonneau, P., Clilverd, M. A., Dudok de Wit, T., Haberreiter, M., Hendry, A., Jackman, C. H., Kretzschmar, M., Kruschke, T., Kunze, M., Langematz, U., Marsh, D. R., Maycock, A. C., Misios, S., Rodger, C. J., Scaife, A. A., Seppälä, A., Shangguan, M., Sinnhuber, M., Tourpali, K., Usoskin, I., van de Kamp, M., Verronen, P. T., and Versick, S.: Solar forcing for CMIP6 (v3.2), Geosci. Model Dev., 10, 2247–2302,, 2017. a

Mauritsen, T., Stevens, B., Roeckner, E., Crueger, T., Esch, M., Giorgetta, M., Haak, H., Jungclaus, J., Klocke, D., Matei, D., Mikolajewicz, U., Notz, D., Pincus, R., Schmidt, H., and Tomassini, L.: Tuning the climate of a global model, J. Adv. Model. Earth Sy., 4, 1–18​​​​​​​,, 2012. a

McManus, J. F., Francois, R., Gherardi, J.-M., Keigwin, L. D., and Brown-Leger, S.: Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes., Nature, 428, 834–837​​​​​​​,, 2004. a

McPhee, M. G.: Turbulent heat flux in the upper ocean under sea ice, J. Geophys. Res., 97, 5365–5379​​​​​​​,, 1992. a, b

Meier, W. N., Fetterer, F., Windnagel, A., and Stewart, J.: NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration, Version 4, National Snow & Ice Data Center [data set],, 2021. a, b

Meinshausen, M., Vogel, E., Nauels, A., Lorbacher, K., Meinshausen, N., Etheridge, D. M., Fraser, P. J., Montzka, S. A., Rayner, P. J., Trudinger, C. M., Krummel, P. B., Beyerle, U., Canadell, J. G., Daniel, J. S., Enting, I. G., Law, R. M., Lunder, C. R., O'Doherty, S., Prinn, R. G., Reimann, S., Rubino, M., Velders, G. J. M., Vollmer, M. K., Wang, R. H. J., and Weiss, R.: Historical greenhouse gas concentrations for climate modelling (CMIP6), Geosci. Model Dev., 10, 2057–2116,, 2017. a

Millero, F. J. and Poisson, A.: International one-atmosphere equation of state of seawater, Deep-Sea Res. Pt. I, 28, 625–629,, 1981. a, b, c

Montoya, M., Griesel, A., Levermann, A., Mignot, J., Hofmann, M., Ganopolski, A., and Rahmstorf, S.: The earth system model of intermediate complexity CLIMBER-3α. Part I: Description and performance for present-day conditions, Clim. Dynam., 25, 237–263,, 2005. a

Morice, C. P., Kennedy, J. J., Rayner, N. A., and Jones, P. D.: Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates: The HadCRUT4 data set, J. Geophys. Res.-Atmos., 117, 1–22​​​​​​​,, 2012. a, b

Müller, S. A., Joos, F., Edwards, N. R., and Stocker, T. F.: Water Mass Distribution and Ventilation Time Scales in a Cost-Efficient, Three-Dimensional Ocean Model, J. Climate, 19, 5479–5499,, 2006. a, b, c, d

Myhre, G., Highwood, E. J., Shine, K. P., and Stordal, F.: New estimates of radiative forcing due to well mixed greenhouse gases, Geophys. Res. Lett., 25, 2715–2718,, 1998. a

Nadeau, L. P., Ferrari, R., and Jansen, M. F.: Antarctic sea ice control on the depth of North Atlantic deep water, J. Climate, 32, 2537–2551,, 2019. a

Nijsse, F. J. M. M., Cox, P. M., and Williamson, M. S.: Emergent constraints on transient climate response (TCR) and equilibrium climate sensitivity (ECS) from historical warming in CMIP5 and CMIP6 models, Earth Syst. Dynam., 11, 737–750,, 2020. a

Niu, G. Y. and Yang, Z. L.: An observation-based formulation of snow cover fraction and its evaluation over large North American river basins, J. Geophys. Res.-Atmos., 112, 1–14​​​​​​​,, 2007. a

Oleson, K. W., Lawrence, D. M., Bonan, G. B., Drewniak, B., Huang, M., Charles, D., Levis, S., Li, F., Riley, W. J., Zachary, M., Swenson, S. C., Thornton, P. E., Bozbiyik, A., Fisher, R., Heald, C. L., Kluzek, E., Lamarque, F., Lawrence, P. J., Leung, L. R., Muszala, S., Ricciuto, D. M., Sacks, W., Sun, Y., Tang, J., and Yang, Z.-L.: Technical Description of version 4.5 of the Community Land Model (CLM) Coordinating, Tech. Rep., No. NCAR/TN-503+STR,, 2013. a

Orr, J. C., Najjar, R. G., Aumont, O., Bopp, L., Bullister, J. L., Danabasoglu, G., Doney, S. C., Dunne, J. P., Dutay, J.-C., Graven, H., Griffies, S. M., John, J. G., Joos, F., Levin, I., Lindsay, K., Matear, R. J., McKinley, G. A., Mouchet, A., Oschlies, A., Romanou, A., Schlitzer, R., Tagliabue, A., Tanhua, T., and Yool, A.: Biogeochemical protocols and diagnostics for the CMIP6 Ocean Model Intercomparison Project (OMIP), Geosci. Model Dev., 10, 2169–2199,, 2017. a

Paul, A., Mulitza, S., Stein, R., and Werner, M.: A global climatology of the ocean surface during the Last Glacial Maximum mapped on a regular grid (GLOMAP), Clim. Past, 17, 805–824,, 2021. a, b

Pedro, J. B., Jochum, M., Buizert, C., He, F., Barker, S., and Rasmussen, S. O.: Beyond the bipolar seesaw: Toward a process understanding of interhemispheric coupling, Quaternary Sci. Rev., 192, 27–46,, 2018. a

Pendergrass, A. G., Conley, A., and Vitt, F. M.: Surface and top-of-atmosphere radiative feedback kernels for CESM-CAM5, Earth Syst. Sci. Data, 10, 317–324,, 2018. a, b, c, d, e

Petoukhov, V., Ganopolski, A., Brovkin, V., Claussen, M., Eliseev, A., Kubatzki, C., and Rahmstorf, S.: CLIMBER-2: a climate system model of intermediate complexity. Part I: model description and performance for present climate, Clim. Dynam., 16, 1–17​​​​​​​,, 2000. a, b, c, d, e, f, g, h, i, j, k

Petoukhov, V., Ganopolski, A., and Claussen, M.: Potsdam – A set of atmosphere statistical-dynamical models: Theoretical background, Tech. Rep., 81, Potsdam Institute for Climate Impact Research, Potsdam, Germany, 2003. a, b, c

Pinardi, N., Rosati, A., and Pacanowski, R. C.: The sea surface pressure formulation of rigid lid models. Implications for altimetric data assimilation studies, J. Marine Syst., 6, 109–119,, 1995. a

Planchon, O. and Darboux, F.: A fast, simple and versatile algorithm to fill the depressions of digital elevation models, Catena, 46, 159–176,, 2002. a

Rahmstorf, S.: A fast and complete convection scheme for ocean models, Ocean Model., 101, 9–11​​​​​​​, 1993. a, b

Rahmstorf, S.: Bifurcations of the Atlantic thermohaline circulation in response to changes in the hydrological cycle, Nature, 378, 145–149,, 1995. a

Rahmstorf, S., Crucifix, M., Ganopolski, A., Goosse, H., Kamenkovich, I., Knutti, R., Lohmann, G., Marsh, R., Mysak, L. A., Wang, Z., and Weaver, A. J.: Thermohaline circulation hysteresis: A model intercomparison, Geophys. Res. Lett., 32, L23605,, 2005. a, b

Redi, M. H.: Oceanic isopycnal mixing by coordinate rotation, J. Phys. Oceanogr., 12, 1154–1158,<1154:OIMBCR>2.0.CO;2, 1982. a

Ritz, S. P., Stocker, T. F., and Joos, F.: A coupled dynamical ocean-energy balance atmosphere model for paleoclimate studies, J. Climate, 24, 349–375,, 2011. a

Robinson, A. and Perrette, M.: NCIO 1.0: a simple Fortran NetCDF interface, Geosci. Model Dev., 8, 1877–1883,, 2015. a

Robinson, A., Alvarez-Solas, J., Montoya, M., Goelzer, H., Greve, R., and Ritz, C.: Description and validation of the ice-sheet model Yelmo (version 1.0), Geosci. Model Dev., 13, 2805–2823,, 2020. a

Roesch, A., Wild, M., Gilgen, H., and Ohmura, A.: A new snow cover fraction parameterization for the ECHAM4 GCM, Clim. Dynam., 17, 933–946,, 2001. a, b

Rossow, W. B. and Schiffer, R. A.: Advances in Understanding Clouds from ISCCP, B. Am. Meteorol. Soc., 80, 2261–2287,<2261:AIUCFI>2.0.CO;2, 1999. a, b

Schaffer, J., Timmermann, R., Arndt, J. E., Kristensen, S. S., Mayer, C., Morlighem, M., and Steinhage, D.: A global, high-resolution data set of ice sheet topography, cavity geometry, and ocean bathymetry, Earth Syst. Sci. Data, 8, 543–557,, 2016. a

Schmidt, G. A., Ruedy, R. A., Miller, R. L., and Lacis, A. A.: Attribution of the present-day total greenhouse effect, J. Geophys. Res.-Atmos., 115, 1–6​​​​​​​,, 2010. a, b, c, d, e

Semtner, A. J.: A Model for the Thermodynamic Growth of Sea Ice in Numerical Investigations of Climate, J. Phys. Oceanogr., 6, 379–389,<0379:AMFTTG>2.0.CO;2, 1976. a, b, c

Shell, K. M., Kiehl, J. T., and Shields, C. A.: Using the radiative kernel technique to calculate climate feedbacks in NCAR's Community Atmospheric Model, J. Climate, 21, 2269–2282,, 2008. a, b, c, d, e

Shin, S. I., Liu, Z., Otto-Bliesner, B. L., Kutzbach, J. E., and Vavrus, S. J.: Southern Ocean sea-ice control of the glacial North Atlantic thermohaline circulation, Geophys. Res. Lett., 30, 68–71,, 2003. a

Smith, C. J., Kramer, R. J., Myhre, G., Forster, P. M., Soden, B. J., Andrews, T., Boucher, O., Faluvegi, G., Fläschner, D., Hodnebrog, Kasoar, M., Kharin, V., Kirkevåg, A., Lamarque, J. F., Mülmenstädt, J., Olivié, D., Richardson, T., Samset, B. H., Shindell, D., Stier, P., Takemura, T., Voulgarakis, A., and Watson-Parris, D.: Understanding Rapid Adjustments to Diverse Forcing Agents, Geophys. Res. Lett., 45, 12023–12031,, 2018. a, b, c, d

Smith, R. S., Gregory, J. M., and Osprey, A.: A description of the FAMOUS (version XDBUA) climate model and control run, Geosci. Model Dev., 1, 53–68,, 2008. a

Stommel, H.: Thermohaline Convection with Two Stable Regimes of Flow, Tellus, 13, 224–230,, 1961. a

Stouffer, R. J. and Manabe, S.: Equilibrium response of thermohaline circulation to large changes in atmospheric CO2 concentration, Clim. Dynam., 20, 759–773,, 2003. a

Subin, Z. M., Riley, W. J., and Mironov, D.: An improved lake model for climate simulations: Model structure, evaluation, and sensitivity analyses in CESM1, J. Adv. Model. Earth Sy., 4, 1–27​​​​​​​,, 2012. a

Tarasov, L., Dyke, A. S., Neal, R. M., and Peltier, W. R.: A data-calibrated distribution of deglacial chronologies for the North American ice complex from glaciological modeling, Earth Planet. Sc. Lett., 315–316, 30–40​​​​​​​,, 2012. a

Tarnocai, C., Canadell, J. G., Schuur, E. a. G., Kuhry, P., Mazhitova, G., and Zimov, S.: Soil organic carbon pools in the northern circumpolar permafrost region, Global Biogeochem. Cy., 23, 2​​​​​​​,, 2009. a

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485–498,, 2012. a

Tierney, J. E., Zhu, J., King, J., Malevich, S. B., Hakim, G. J., and Poulsen, C. J.: Glacial cooling and climate sensitivity revisited, Nature, 584, 569–573,, 2020. a, b, c

Trenberth, K. E. and Caron, J. M.: Estimates of Meridional Atmosphere and Ocean Heat Transports, J. Climate, 14, 3433–3443,<3433:EOMAAO>2.0.CO;2, 2001. a

Trenberth, K. E., Smith, L., Qian, T., Dai, A., and Fasullo, J.: Estimates of the Global Water Budget and Its Annual Cycle Using Observational and Model Data, J. Hydrometeorol., 8, 758–769,, 2007. a

Vavrus, S. and Waliser, D.: An improved parameterization for simulating Arctic cloud amount in the CCSM3 climate model, J. Climate, 21, 5673–5687,, 2008. a

Vellinga, M. and Wood, R. A.: Global climatic impacts of a collapse of the atlantic thermohaline circulation, Climatic Change, 54, 251–267,, 2002. a

Weaver, A. J., Eby, M., Wiebe, E. C., Ewen, T. L., Fanning, A. F., MacFadyen, A., Matthews, H. D., Meissner, K. J., Saenko, O., Schmittner, A., Yoshimori, M., Bitz, C. M., Holland, M. M., Duffy, P. B., and Wang, H.: The UVic earth system climate model: Model description, climatology, and applications to past, present and future climates, Atmos. Ocean, 39, 361–428,, 2001. a, b, c

Weber, S. L., Drijfhout, S. S., Abe-Ouchi, A., Crucifix, M., Eby, M., Ganopolski, A., Murakami, S., Otto-Bliesner, B., and Peltier, W. R.: The modern and glacial overturning circulation in the Atlantic ocean in PMIP coupled model simulations, Clim. Past, 3, 51–64,, 2007. a

Weijer, W., Cheng, W., Drijfhout, S. S., Fedorov, A. V., Hu, A., Jackson, L. C., Liu, W., McDonagh, E. L., Mecking, J. V., and Zhang, J.: Stability of the Atlantic Meridional Overturning Circulation: A Review and Synthesis, J. Geophys. Res.-Oceans, 124, 5336–5375,, 2019. a

Wetherald, R. T. and Manabe, S.: Cloud Feedback Processes in a General Circulation Model, J. Atmos. Sci., 45, 1397–1416,<1397:CFPIAG>2.0.CO;2, 1988. a

Wild, M., Folini, D., Schär, C., Loeb, N., Dutton, E. G., and König-Langlo, G.: The global energy balance from a surface perspective, Clim. Dynam., 40, 3107–3134,, 2013. a

Willeit, M.: CLIMBER-X v1.0, Zenodo [code],, 2022. a

Willeit, M. and Ganopolski, A.: PALADYN v1.0, a comprehensive land surface–vegetation–carbon cycle model of intermediate complexity, Geosci. Model Dev., 9, 3817–3857,, 2016. a, b, c, d, e, f, g

Willeit, M., Ganopolski, A., and Feulner, G.: Asymmetry and uncertainties in biogeophysical climate–vegetation feedback over a range of CO2 forcings, Biogeosciences, 11, 17–32,, 2014. a

Willeit, M., Ganopolski, A., Calov, R., and Brovkin, V.: Mid-Pleistocene transition in glacial cycles explained by declining CO2 and regolith removal, Science Advances, 5, eaav7337,, 2019. a

Yamamoto, G. and Tanaka, M.: Increase of Global Albedo Due to Air Pollution, J. Atmos. Sci., 29, 1405–1412,<1405:IOGADT>2.0.CO;2, 1972. a

Yang, H., Li, Q., Wang, K., Sun, Y., and Sun, D.: Decomposing the meridional heat transport in the climate system, Clim. Dynam., 44, 2751–2768,, 2015. a

Yin, J., Stouffer, R. J., Spelman, M. J., and Griffies, S. M.: Evaluating the uncertainty induced by the virtual salt flux assumption in climate simulations and future projections, J. Climate, 23, 80–96,, 2010. a

Zalesak, S. T.: Fully multidimensional flux-corrected transport algorithms for fluids, J. Comput. Phys., 31, 335–362,, 1979. a, b, c

Zelinka, M. D., Klein, S. A., and Hartmann, D. L.: Computing and partitioning cloud feedbacks using cloud property histograms. Part II: Attribution to changes in cloud amount, altitude, and optical depth, J. Climate, 25, 3736–3754,, 2012.  a

Zelinka, M. D., Myers, T. A., McCoy, D. T., Po-Chedley, S., Caldwell, P. M., Ceppi, P., Klein, S. A., and Taylor, K. E.: Causes of Higher Climate Sensitivity in CMIP6 Models, Geophys. Res. Lett., 47, e2019GL085782​​​​​​​,, 2020. a

Zika, J. D., Skliris, N., Blaker, A. T., Marsh, R., Nurser, A. J., and Josey, S. A.: Improved estimates of water cycle change from ocean salinity: The key role of ocean warming, Environ. Res. Lett., 13, 074036,, 2018. a, b

Short summary
In this paper we present the climate component of the newly developed fast Earth system model CLIMBER-X. It has a horizontal resolution of 5°x5° and is designed to simulate the evolution of the Earth system on temporal scales ranging from decades to >100 000 years. CLIMBER-X is available as open-source code and is expected to be a useful tool for studying past climate changes and for the investigation of the long-term future evolution of the climate.