Evaluation and climate sensitivity of the PlaSim v.17 Earth System Model coupled with ocean model components of different complexity

A set of experiments is performed with coupled atmosphere-ocean configurations of the Planet Simulator, an Earth-system Model of Intermediate Complexity (EMIC), in order to identify under which set of parameters the model output better agrees with observations and reanalyses of the present climate. Different model configurations are explored, in which the atmospheric module of PlaSim is coupled with two possible ocean models, either a simple mixed-layer (ML) ocean with a diffusive trans5 port parameterization or a more complex dynamical Large-Scale Geostrophic (LSG) ocean, together with a sea-ice module. In order to achieve a more realistic representation of present-day climate, we performed a preliminary tuning of the oceanic horizontal diffusion coefficient for the ML ocean and of the vertical oceanic diffusion profile when using LSG. Model runs under present-day conditions are compared, in terms of surface air temperature, sea surface temperature, sea ice cover, precipitation, radiation fluxes, ocean circulation, with a reference climate from observations and reanalyses. Our results indicate 10 that, in all configurations, coupled PlaSim configurations are able to reproduce the main characteristics of the climate system, with the exception of the Southern Ocean region in the PlaSim-LSG model, where surface air and sea surface temperatures are warm-biased and sea ice cover is by consequence highly underestimated. The resulting sets of tuned parameters are used to perform a series of model equilibrium climate sensitivity (ECS) experiments, with the aim to identify the main mechanisms contributing to differences between the different configurations and 15 leading to elevated values of ECS. In fact, high resulting global ECS values are found, positioned in the upper range of CMIP5 and recent CMIP6 estimates. Our analysis shows that a significant contribution to ECS is given by the sea-ice feedback mechanisms and by details of the parameterization of meridional oceanic heat transport. In particular, the configurations using a diffusive heat transport in the mixed layer present an important sensitivity in terms of radiative forcing to changes in sea-ice cover, leading to an important contribution of sea-ice feedback mechanisms to ECS. 20 1 https://doi.org/10.5194/gmd-2020-245 Preprint. Discussion started: 28 October 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
While numerical Global Climate Models (GCMs) are continuously improved to understand the main processes at work in the climate system and to evaluate the impacts of changes in radiative forcing, the complexity of these models has considerably grown over the past decades and simpler models have remained attractive for addressing research questions, understanding specific processes at work and for academic education. In these models atmospheric and climate processes are simplified or highly parameterized so that key mechanisms become easier to study. In particular, Earth-system Models of Intermediate al. (2020); Meehl et al. (2020)). Two processes which we find to have a significant impact on equilibrium climate sensitivity, sea-ice feedbacks and meridional oceanic heat transport, are investigated in more detail.
The paper is structured as follows: in Section 2 the various components of the PlaSim model, the characteristics of the ML and LSG ocean models and their coupling with the atmospheric component are described; Section 3 presents the tuning of ocean parameters in the PlaSim-ML and PlaSim-LSG configurations; in Section 4 the simulated climate is evaluated, including an assessment of the energy and water balance of the model; Section 5 deals with the equilibrium climate sensitivity of the 5 model in the selected configurations and Section 6 summarizes the main results and concludes the paper. In this section, the main characteristics of the PlaSim Earth system Model of Intermediate Complexity  are briefly illustrated, but we refer to the PlaSim reference manual, Lunkeit et al. (2011), for further details. The dynamical core 10 of PlaSim is a simplified General Circulation Model, the Portable University Model of Atmosphere (PUMA), based on the moist primitive equations representing the conservation of momentum, mass and energy (Fraedrich et al., 2005a) and using spectral methods to numerically solve them (Orszag, 1970;Eliasen et al., 1970). In the vertical, a σ-coordinate system and a finite-difference method to solve equations are used. The equations are time integrated with a leap-frog semi-implicit time stepping scheme with time filter (Hoskins & Simmons, 1975;Simmons et al., 1978;Robert, 1981;Asselin, 1972).
suitable for climate studies under conditions far from present day, we will not focus on this configuration in this study. Instead oceanic transports are parameterized by the addition of an horizontal diffusion term to the temperature equation: where T M L is the mixed-layer temperature, F a describes the net energy exchanges with the atmosphere and K h is a horizontal temperature diffusion coefficient (with a starting default value of 1000 m 2 s −1 , but we tuned this value in our experiments).
The sea-ice distribution can either be prescribed by climatology or simulated by a thermodynamic sea-ice model based on 5 the zero layer model by Semtner (1976), which computes the thickness of sea ice from the thermodynamic balance at the top and at the bottom of the sea-ice layer. This model assumes a linear temperature gradient in the ice and prevents ice from storing heat. Sea ice is formed if the ocean temperature drops below the freezing point (set to 271.25 K) and is melted if the ocean temperature exceeds that value .

The Large Scale Geostrophic ocean circulation model 10
The Large Scale Geostrophic (LSG) ocean circulation model (Maier-Reimer et al., 1993;Drijfhout et al., 1996) is based on the primitive equations in a three-dimensional system, including the momentum equation, the continuity equation describing conservation of water and salinity, the thermodynamic equation with salinity. Please see the Large Scale Geostrophic Model report, Maier-Reimer & Mikolajewicz (1992), for details. The model is based on the observation that, since for a large scale ocean circulation model developed for climate studies the characteristic spatial scales are large compared with the internal 15 Rossby radius of deformation and the characteristic temporal scales are large compared with the periods of gravity modes and barotropic Rossby wave modes (Hasselmann, 1982), the nonlinear terms in the Navier-Stokes equations can be neglected.
Furthermore the vertical friction is neglected and the hydrostatic and the Boussinesq approximations are applied (Maier-Reimer & Mikolajewicz, 1992).
Turbulent motions are parameterized by means of a vertical oceanic diffusion coefficient, A v , which is a rather simple 20 function of the vertical coordinate, z (Bryan & Lewis, 1979): where a * is the vertical diffusion coefficient at a reference depth z * , a range defines the considered depth range from the surface to the bottom, and λ is the rate at which the vertical diffusion coefficient varies with depth near z * .
A long time step of 10 days is permitted by the implicit time integration scheme. The model has two staggered 5°x 5°25 horizontal grids (yelding an effective grid resolution of 3.5°), so that the variables of the model are defined on a semi-staggered E-type grid (Arakawa & Lamb, 1977). The components of horizontal velocity and the wind stress are defined on "vector points", while potential temperature, salinity, heat and freshwater fluxes, sea-surface height, pressure and vertical velocity are defined on "scalar points". The depth of the scalar points is usually defined as the maximum depth of the four surrounding vector points. interpolation ensures global conservation of energy and water Lorenz, 2006 A full ocean step is performed in which LSG calculates the ocean heat flux due to advective (advection, horizontal diffusion) and convective (vertical transport, vertical diffusion, convective adjustments) processes. The resulting distribution of the upper layer temperature after ∆t (LSG) determines the ocean heat flux, which is then given back to PlaSim. Equation. 3 is further 15 modified to take into account small differences in the ML depth used by PlaSim and the upper layer of LSG and to correctly close the energy balance (see Lorenz (2006) for details).
Sea ice in LSG is prescribed as calculated in the ML module of PlaSim. When calculating over hundreds of years, sea ice grows unconstrained in some isolated grid-points in the Antarctic ocean. In order to avoid the consequent increase of salinity of open water in the upper layer of LSG, sea-ice thickness is forced not to exceed 9 m.

20
Furthermore, the atmospheric wind stress and the freshwater flux (with a constant annual mean flux correction) are averaged over the coupling interval before they are transferred to the ocean.

Coupled ocean tuning
In the following we focus on the tuning of two oceanic parameters which we have found to be crucial in determining the ability of the model in reproducing the observed reference climate. In particular we consider: i) the horizontal diffusion coefficient 25 (when using the ML ocean) and ii) the vertical diffusion coefficient (when using the LSG ocean).

Mixed-layer ocean horizontal diffusion
The simulations with PlaSim-ML were performed including the dynamic sea-ice module, at T21 spatial resolution and with a fixed atmospheric CO 2 concentration (perpetual simulations) set to 354 ppm (representing conditions in 1990). Please notice that in the following we are comparing results from simulations with fixed CO 2 concentrations with observed climatologies in the years 2005-2015: using CO 2 concentrations from a previous period (1990) aims at compensating climate warming experienced by a climate model under fixed greenhouse gas forcing, assuming a feedback parameter of about 1 Wm −2 K −1 5 and the current radiative imbalance at TOA of about 0.5 Wm −2 . We performed six 60 year long runs, with fixed atmospheric CO 2 concentrations, each having a different horizontal diffusion coefficient, K h , which was varied in the range from 10 3 m 2 s −1 to 10 6 m 2 s −1 . Panel a of Fig. 1 shows the global mean of surface air temperature time series for all these simulations, each colored line representing one simulation performed with a different value of K h .  (Dee et al., 2011), to be compared with the model outputs.
In order to identify which K h value(s) leads to a better agreement between the model and the observations, we first computed, for each perpetual simulation, a time average over the last 30 years and then we compared the zonally-averaged surface air temperature in the model and in the ERA-Interim dataset, as shown in Fig. 1 panel b. The figure clearly suggests that by 5 choosing a single value of K h for the entire globe leads to either a cold bias in the Northern Hemisphere (NH) or a warm bias in the Southern Hemisphere (SH). In particular, the zonal mean temperatures in the simulation performed with K h = 10 5 m 2 s −1 (green line) are in a good agreement with the reanalysis data in the NH, but give rise to a warm bias (up to 21 K) in the SH.
Using instead K h = 10 4 m 2 s −1 (red line), simulated surface air temperatures are in good agreement with the reanalyses in the SH but they are too cold (up to 20 K) in the NH. This suggests the use of two different diffusion coefficients, to be 10 separately applied in the Northern and Southern Hemisphere, i.e. K h = 10 5 m 2 s −1 in the NH and K h = 10 4 m 2 s −1 in the SH (this choice can be, to some extent, justified physically by the observed differences in meridional heat transport between the two hemispheres, particularly the strong North-South asymmetry observed in the Atlantic basin). Panel c of Fig. 1 shows the results of this simulation, giving rise to a simulated zonally-averaged surface air temperature in very good agreement with the reanalyses. The maximum difference between the model results and the ERA-Interim values is about 3 K. 15 We performed the same experiments (PlaSim-ML with dynamic sea ice, CO 2 = 354 ppm) running the model at higher resolution, T42. In this case, the horizontal diffusion coefficients K h has been tuned testing five different values, ranging from 10 4 m 2 s −1 to 10 5 m 2 s −1 . In Fig. 1, panel d and e show the results of these experiments in terms of, respectively, surface air temperature time series and zonally-averaged surface air temperature for each simulation performed with a different K h value.
Also in this case (T42 resolution) a single diffusion coefficient does not allow to find a good agreement between the model 20 and the reanalyses at the global level. In particular, panel e would suggest to use K h = 10 5 m 2 s −1 (green line) in the NH and K h = 3 · 10 4 m 2 s −1 (red line) in the SH. The results of the simulation performed with the two coefficients are shown in panel f .

Large Scale Geostrophic ocean vertical diffusion
We performed the series of experiments with PlaSim coupled with the Large Scale Geostrophic ocean model using dynamic 25 sea ice, a fixed atmospheric CO 2 concentration equal to 354 ppm and atmospheric resolution T21. We initially performed two perpetual, 2000 year long runs, each having a different vertical diffusion profile. For the first run (run 1), we used the vertical diffusion coefficient profile suggested by Bryan & Lewis (1979) Table 1 and Fig. 2, blue line). For the second run (run 2), we used modified parameters which were found by Sciascia (2008) in ocean-only tuning experiments, with a surface value of A v = 0.8 · 10 −4 m 2 s −1 30 and a bottom value of A v = 1.3 · 10 −4 m 2 s −1 (Table 1 and Fig. 2, magenta line). Table 1. Vertical diffusion parameters used in LSG model simulations. Av n. 1, 2 and 3 correspond respectively to the blue, magenta and green profiles shown in Fig. 2.
run 1 0.7958 · 10 −4 0.3345 · 10 −4 2500 4.5 · 10 −3 run 2 1.0479 · 10 −4 0.1673 · 10 −4 2500 4.5 · 10 −3 run 3 0.8714 · 10 −4 0.2843 · 10 −4 2500 4.5 · 10 −3 using these two coefficients. We computed the time mean over the last 1000 years and we analysed the zonal mean of surface air temperature (panel c) and the anomaly with respect to the ERA Interim dataset (inset box). In the Northern Hemisphere, the simulated temperature is negatively biased with coefficient n. 1 due to the AMOC collapse, while using coefficient n. 2 the PlaSim-LSG model maintains the oceanic circulation but overestimates (up to 12 K) surface air temperatures from 40°to 90°i n the Southern Hemisphere. plays an important role, while its variation below 2000 m has no significant impact. Based on these results, we chose a vertical diffusion profile (run 3) with a surface value of A v = 0.45 · 10 −4 m 2 s −1 and a bottom value of A v = 1.3 · 10 −4 m 2 s −1 (see Table 1 and Fig. 2, green line), which is the best compromise between active AMOC (about 19 Sv) and lower temperature in the Southern Hemisphere (see Fig. 3). In this configuration, the coupled model PlaSim-LSG reaches an equilibrium state after about 500 years of simulation and the zonally-averaged maximum warm bias in the Southern Ocean is 8 K. In the attempt to reduce this bias, we have modified other parameters (within physically acceptable limits), such as the cloud albedo, the 5 continental ice-sheet albedo, the oceanic albedo and horizontal diffusion, the ozone concentration in the atmosphere, but with a negligible improvement of the resulting climate. Based on these results, we concluded that the simulation with profile n. 3 (green profile in Fig. 2) best reproduces temperature estimates from reanalyses and we chose it for the following PlaSim-LSG reference run.

Model climate 10
The preliminary tuning of ocean parameters described in the previous section allowed us to define three reference configurations which in the following we use to perform equilibrium climate sensitivity calculations. The first configuration consists of PlaSim coupled with the mixed-layer ocean, run at T21 spatial resolution and with two different horizontal oceanic diffusion coefficients, for the Northern Hemisphere (K h = 10 5 m 2 s −1 ) and for the Southern Hemisphere (K h = 10 4 m 2 s −1 ). The second configuration is similar to the first one but the spatial resolution is finer (T42) and with a different horizontal diffusion oceanic   PlaSim-LSG (panel n) better correlates with the observed precipitation than PlaSim-ML. Finally, TOA net fluxes (panels o-q) show a negative anomaly in the tropics and subtropics, where the upward radiation is overestimated, and a positive anomaly in the mid-latitude and polar zones, where the upward radiation is underestimated. Furthermore, a positive anomaly of net 20 radiation is observed on the western coast of America and Africa, consistently with the warm biased surface air temperatures.

Energy balance
For any coupled climate model to reach a stationary stable state, the TOA and the surface energy balances should be close to zero (neglecting geothermal heating) since the model should present no significant internal energy sources or sinks. Table   3 compares the global energy balance resulting from our reference simulations with estimates from Stephens et al. (2012).
While the former are the results of perennial model simulations with PlaSim, the latter refers to an observed climatology . We tested if this imbalance is caused by a missing conservation of water mass in the model atmosphere (possibly to transport errors), but the absolute value of the global average freshwater flux P-E is smaller than 10 −4 mm/day, equivalent to a very small latent heat flux smaller than 2.6 · 10 −3 Wm −2 for all tested PlaSim configurations, indicating that water is well conserved in the PlaSim atmosphere. Also, all ML simulations present a negative net energy flux at the surface, suggesting some non-conservation of energy (equivalent to an energy production) in the mixed-layer ocean. Both the TOA- 15 surface imbalance and the net surface flux bias are reduced in the T42 ML simulation, suggesting that these biases may be resolution-dependent. Overall these energy imbalances are small compared to those reported for CMIP5 and CMIP3 models (which could exceed 1 W/m 2 in magnitude at TOA (Mauritsen et al., 2012)). Equilibrium climate sensitivity (ECS) is defined as the equilibrium change in global mean surface air temperature after an instantaneous doubling of atmospheric CO 2 relative to pre-industrial levels (IPCC, 2013). Climate sensitivity can be diagnosed following the approach by Gregory et al. (2004) and here we apply this method to PlaSim-ML and PlaSim-LSG simulations.
When a radiative forcing R (Wm −2 ) is applied to the model, the model responds with a change in the net TOA radiative flux ∆F (Wm −2 ) and, in order to restore the radiative equilibrium, the global mean surface air temperature, ∆T , changes, until 5 ∆F is returned to zero. R, ∆F and ∆T are related by the following equation: where λ (Wm −2 K −1 ) is referred to as climate feedback parameter. If ∆F is assumed to be a linear function of ∆T , both the radiative forcing and the feedback parameter can be diagnosed by linear regression: R is the intercept at ∆T = 0 and λ is the slope (multiplied by -1). The equilibrium temperature change can be estimated extrapolating the heat balance to equilibrium, 10 that is ∆F = 0 and ∆T eq = R/λ. If the forcing is a doubling of CO 2 , ∆T eq is the equilibrium climate sensitivity by definition.
We performed a first set of simulations using dynamic sea ice (subscript d in subsequent text): the first part of each simulation is a perennial run with pre-industrial boundary conditions, so the CO 2 concentration in the atmosphere is set to 285 ppm (1xCO 2 ); the second part is a perennial run in which the CO 2 concentration is instantaneously increased at 1.5, 2, 3 or 4 times the value of the pre-industrial simulation. These simulations were made with the three reference configurations of PlaSim.

15
Each half-simulation is 100 years long when using PlaSim-ML and 2000 years long when using PlaSim-LSG. The yellow and the red lines in Fig. 5 show the change in net TOA radiative flux versus the change in global mean surface air temperature for each dynamic-ice simulation with doubled (2xCO 2 ) and quadrupled (4xCO 2 ) atmospheric CO 2 . Changes are computed with respect to the corresponding 1xCO 2 part of the simulation.  We derive the estimates of R d (intercept) and λ d (slope multiplied by -1) through ordinary least squares regression and 20 we compute ∆T eq d = R d /λ d for all the PlaSim configurations (Table 4). The confidence intervals are obtained as standard deviation on the parameter estimates and error propagation. The resulting ECS for dynamic sea ice is 6.23 K using PlaSim-ML T21, 5.45 K using PlaSim-ML T42 and 4.26 K using PlaSim-LSG T21, using the results from the CO 2 doubling experiments.
In Fig. 6 we compare these results with values from other models. In particular, the grey boxplots give an indication of the distribution of CMIP5 values (the whiskers extend to the highest and lowest data) discussed in Andrews et al. (2012). Radiative 5 forcing and climate feedback values of our model are within the range estimated for CMIP5 models, but only the PlaSim-LSG coupled model gives an equilibrium climate sensitivity within the CMIP5 range (2.1-4.7 K), though close to the upper limit.
The orange boxplot represents the ECS values found in other EMICs (Pfister & Stocker, 2017), which are in good agreement with CMIP5 models but have a wider range of values (1.5-5.5 K). Finally, the blue boxplot shows the most recent range of ECS for CMIP6 models (1.8-5.6 K; Zelinka et al. (2020)). With respect to CMIP5, some components of CMIP6 models have been 10 improved, for example low clouds and shallow convection are better represented (Voldoire et al., 2019) or a more advanced treatment of aerosol is included (Wyser et al., 2020), and stronger positive cloud feedbacks from decreasing extratropical low cloud coverage and albedo have contributed to increased ECS in some of them (Zelinka et al., 2020;Meehl et al., 2020). Since PlaSim does not include such a level of accuracy (for example, it has no parameterization for aerosol-cloud feedbacks), its high climate sensitivity cannot be related to these processes. Our results can also be compared also with an ECS estimate for 15 a modified PlaSim ML configuration at T21 in Ragone et al. (2016), where the very high value of 8.1 K was reported and attributed to the removal of the diurnal and seasonal cycles in the model. Please notice also that, while EMIC ECS values reported in Fig. 6 were obtained, like for PlaSim, from CO 2 doubling experiments, the reported CMIP5 and CMIP6 results were obtained dividing by two the results from 4xCO 2 experiments. As shown in Table 4 for PlaSim and as also reported in Pfister & Stocker (2017) for other EMICS, the ECS values obtained from quadrupling experiments may be lower than those obtained from doubling experiments, although often used without distinction as estimates of ECS in the literature.  (Pfister & Stocker, 2017) and CMIP6 models (Zelinka et al., 2020) . The difference between the two PlaSim-ML configurations and the PlaSim-LSG configuration can partly be explained by 5 features of the ocean circulation. Using PlaSim-LSG, the AMOC is active (about 20 Sv) in the pre-industrial and 1.5xCO 2 runs, but it collapses (less than 5 Sv) in the 2xCO 2 , 3xCO 2 and 4xCO 2 simulations. As shown in Fig. 3b, global average temperatures are affected significantly by the state of AMOC, with a cooling of up to 1°C in runs with a shutdown of AMOC.
As a consequence, the equilibrium climate sensitivity in such runs is smaller than it would be if the AMOC had remained active. Unlike PlaSim-LSG and CMIP models, the PlaSim-ML configurations doesn't include an AMOC representation, so it 10 can't weaken and this could contribute to the reported higher ECS.
The relatively high values of ECS found for PlaSim are related to low values of the feedback parameter λ d . We determined that an important contribution can be traced also to elevated values in magnitude of the ice-feedback parameter, as we assessed following the approach of Caldeira & Cvijanovic (2014). To this end, we performed a second set of simulations, similar to the first one but with prescribed sea ice (subscript p): twelve climatological monthly ice extents were derived from the pre-15 industrial dynamic-ice simulation and were prescribed in the model. The cyan and the blue lines in Fig. 5 show the change in net TOA radiative flux versus the change in global mean surface air temperature for each prescribed-ice simulation with doubled and quadrupled CO 2 concentration. Using prescribed sea ice, the change in TOA radiative flux at equilibrium is not zero (see Fig. 5) because some energy has to be removed from or added to the system in order to maintain the climatological sea-ice thickness (Caldeira & Cvijanovic, 2014). Also in this case, we computed λ p (slope multiplied by -1) through ordinary linear least squares regression for all the PlaSim configurations (Table 4). The sea-ice feedback parameter λ i (last column of Table 4) is negative in our sign convention and is obtained subtracting the feedback parameter of dynamic-ice simulations from that of prescribed-ice simulations. These results can be compared with the slab-ocean experiments performed by Caldeira & Cvijanovic (2014) with the National Center for Atmospheric Research's Community Earth System Model (CESM): they 5 report a λ i of -0.21 ± 0.19 Wm −2 K −1 in the doubling CO 2 experiments and -0.30 ± 0.06 Wm −2 K −1 in the quadrupling CO 2 experiments. The feedback parameter of sea ice in the PlaSim-ML configurations is significantly higher in absolute value than in CESM, suggesting that sea ice plays an important role in determining the ECS of the model. However, this contribution also depends on the extent of sea ice either in the pre-industrial climate (see I P I d in Table 4) or in the future climates (see I d ). Indeed, the sea-ice area is very different in the three configurations of PlaSim. For example, in the PlaSim-ML model at 10 T42 the pre-industrial sea-ice area is less extended than in the PlaSim-ML model at T21, so the sea-ice contribution to ECS is smaller. Furthermore, we recall that sea ice is almost completely absent in the Southern Ocean using PlaSim-LSG, therefore the sea-ice contribution to ECS is reduced compared to the configurations with the ML ocean. A factor which also contributes to the low values of λ i in the PlaSim-LSG configuration, is the fact that in our simulations in the pre-industrial simulation with prescribed sea ice, the AMOC in the LSG ocean model collapses to very low values after about 1000 years, while using 15 the dynamic sea-ice treatment the AMOC is strong in the pre-industrial climate. The AMOC collapse in the prescribed-ice simulation makes the pre-industrial global temperature colder (as reported above), the slope λ p smaller and the difference λ i closer to zero than it would be with an active AMOC. Therefore in the PlaSim-LSG model the feedback parameter λ i , obtained subtracting dynamic-ice from prescribed-ice simulations, includes not only the sea-ice contribution but also the AMOC effect, which is positive in our sign convention.

20
The high impact of sea-ice related feedbacks in the model cannot be linked to a too high sea-ice albedo in PlaSim: in fact we compared the average sea-ice albedo of PlaSim and that of EC-Earth, a global state-of-the-art climate model with higher complexity and spatial resolution (Hazeleger et al., 2011;Döscher et al., 2020), which has an ECS of 4.3 K in the newer model version (Wyser et al., 2020). The average sea-ice albedo of PlaSim-ML T21 is 0.58 for pixels with more than 99% area coverage, lower than the average sea-ice albedo of EC-Earth (0.80). Since a smaller sea-ice albedo weakens the 25 ice-albedo feedback, we can conclude that the strong impact of dynamic sea ice in PlaSim is not likely due to the ice albedo parameterization employed in the model. In fact, we performed a series of climate sensitivity runs with PlaSim (not shown) in which we modified maximum sea-ice albedo and its dependence on temperature, without finding significantly lower values of the ECS.
A confusing factor between the different model configurations is the fact that they are all characterized by different average 30 sea-ice extents in the starting preindustrial experiments. To better compare the impact of the sea-ice feedback, we can use the same approach as Caldeira & Cvijanovic (2014), comparing the dynamic-ice and prescribed-ice simulations, to define a measure of the radiative forcing associated with changes in sea-ice area. An equivalent formulation with our symbols is: which represents a measure of the radiative forcing that should be provided to prescribed-ice experiments in order to undergo the same global mean surface air temperature change as in dynamic-ice experiments. Figure 7 shows the sea-ice radiative forcing versus the relative sea-ice area (computed as the difference between I p and I d , see Table 4  . Sea-ice radiative forcing (see Eq. 5) as a function of the relative sea-ice area (prescribed-ice area minus dynamic-ice area, see Table 4). Values are obtained by estimating the radiative forcing that should be provided to prescribed-ice experiments in order to undergo the same global mean surface air temperature change as in dynamic-ice experiments.
The configurations of PlaSim with the mixed layer, in which the oceanic transport is parameterized by a horizontal diffusion term, show a radiative forcing associated with changes in sea-ice area which is significantly higher than CESM, consistently with the higher value of λ i in absolute value which we reported above. One important difference between our experimental 10 results and the CESM results is that the latter experiments used a mixed-layer ocean with a prescribed ocean heat flux (Q-flux).
In order to verify the role played by the specific parameterization of oceanic transport, we performed also a series of CO 2 increase experiments with dynamical and prescribed ice, using PlaSim with a Q-flux correction instead of horizontal diffusion in the mixed layer (the Q-flux was derived by the model from present-day AMIP experiment forced with observed SSTs). As shown in the figure, the configuration of PlaSim-ML T21 with the flux correction gives similar results to CESM in terms of ice and as a consequence the ECS of the model. To further explore this point, Fig. 7 also shows three additional simulations which use a single horizontal diffusion coefficient for both the hemispheres, K h = 10 5 m 2 s −1 (yellow line), K h = 10 4 m 2 s −1 (cyan line) and K h = 10 3 m 2 s −1 (magenta line). The first and the second coefficient are those used for the Northern and the Southern Hemispheres in the PlaSim-ML (T21) configuration. We can notice that PlaSim-LSG and PlaSim-ML with the lowest diffusion coefficient (K h = 10 3 m 2 s −1 ) have a slope similar to CESM, while using a higher ocean diffusion, the effective changes in radiative forcing associated with changes in sea-ice area increase for increasing horizontal oceanic diffusion. Therefore 5 ultimately the equilibrium climate sensitivity of the PlaSim-ML model using a diffusive term to represent heat transport, depends crucially also on the choice of K h .

Conclusions
In this paper the evaluation and the climate sensitivity of coupled configurations of the PlaSim model, using two different ocean modules (a simple mixed layer or the Large Scale Geostrophic ocean circulation model), are presented. In order to parameterize heat transport in the ML, a horizontal diffusion coefficient K h was used, a configuration which potentially could be preferable for climate change or paleoclimatic studies, due to its ability to adapt heat transport to varying pole-equator gradients, compared to a fixed Q-flux.
A preliminary tuning to achieve realistic zonally-averaged near-surface temperature profiles suggests the use of two different values for K h , one in the Northern Hemisphere (NH) and another one in the Southern Hemisphere (SH). To the same end, for 15 the coupled PlaSim-LSG model we found that, when using the arctan-shaped vertical ocean diffusion profile following Bryan & Lewis (1979), values ranging from 0.45 · 10 −4 m 2 s −1 at the top to 1.3 · 10 −4 m 2 s −1 at the bottom of the ocean should best be used. Most climatic variables are well simulated and their anomalies are reduced using the finer T42 resolution. There is a warm bias in the Southern Hemisphere (between 40 and 90°S) and the Antarctic sea ice is underestimated in the coupled PlaSim-LSG model. PlaSim, in all tested configurations, does not provide a perfect energy balance, probably due to its relatively 20 coarse spatial resolution. In fact, we found smaller imbalance in the simulations at higher resolution. Still, overall these energy imbalances are small compared to those reported in literature for several other global climate models.
The tuned PlaSim configurations were used to assess the equilibrium climate sensitivity (ECS) of the model from CO 2 doubling experiments with dynamic sea ice and to compare it with values found in CMIP5 models, other state-of-the-art EMICs and CMIP6 models. The ECS of the model is found to be particularly high in the ML configurations, with an ECS 25 of 6.23 K using PlaSim-ML at T21 and 5.45 K using PlaSim-ML at T42, compared to 4.26 K using PlaSim-LSG at T21.
Only the latter value is within the CMIP5 range, though close to the upper limit. One important factor contributing to higher ECS in PlaSim-ML compared to PlaSim-LSG is that unlike the PlaSim-LSG configuration and complex CMIP models with a 3D ocean, in PlaSim with the ML ocean the AMOC can't weaken or collapse and we have seen also in our simulations that reduced AMOC strength will have a cooling effect on the global average temperatures. The high ECS of the model in the ML 30 configurations has also been found to be related to elevated values in magnitude of the ice-feedback parameter, as computed performing prescribed sea-ice simulations. In fact, overall, our climate sensitivity experiments with the PlaSim EMIC reveal that details of the oceanic heat transport play a very important role in determining the sea-ice feedback parameter and as a consequence the ECS of the model. When using a diffusive term in the ML, with values of the horizontal diffusion parameter K h which allow for a realistic meridional temperature distribution in present-day experiments, changes in average sea-ice area have a much stronger radiative impact, compared to very low values of the diffusion coefficient or to using a Q-flux approach to represent transport. Since using a diffusive term may be preferable to a fixed Q-flux in some cases for studying climate responses far from present-day conditions (such as paleoclimatic or eexoplanetary studies), this impact may have to be taken 5 carefully into account. The fact that the PlaSim configuration with LSG presents a lower radiative impact of sea ice, comparable to that of experiments with a Q-flux, reveals that indeed sea-ice feedbacks may be overestimated in the configurations using a ML with a diffusive transport parameterization. An additional factor contributing to varying values of ECS between different configurations is that, due to the relevant impact of sea-ice feedbacks, the sea-ice extent in preindustrial simulations plays an important role, and this contributes to the lower ECS of the configuration using LSG, which is affected by a strong southern-10 ocean bias and almost no sea ice in the southern hemisphere.  Hazeleger, W., Wang, X., Severijns, C.,Ştefanescu, S., Bintanja, R., Sterl, A., Wyser, K., Semmler, T., Yang, S., Van den Hurk, B., Van Noije, T., Van der Linden, E., Van der Wiel, K.: EC-Earth V2.2: description and validation of a new seamless earth system prediction model,