Abstract Equilibrium climate sensitivity (ECS) – the global temperature response to doubling CO 2 – can be better understood by examining past climate states and the feedback mechanisms regulating CO 2 -induced warming. While studies agree that ECS increases with higher CO 2 , the range of Earth’s historical ECS and its underlying drivers remain incompletely understood. Here, we use slab ocean Community Earth System Model simulations to analyze four distinct periods in Earth's climate history with substantially different continental configurations: the late Cretaceous, early Eocene, late Oligocene, and preindustrial. Our results show that ECS varies by over 2.5°C, ranging from 4.04°C to 6.66°C. We analyze the contributions of CO 2 background, geography, and ocean heat transport to ECS variability and decompose the total climate feedback parameter into the water vapor, cloud, surface albedo, and temperature feedbacks. Using a consistent model framework across geological epochs, we provide new constraints on ECS sensitivity to boundary conditions and offer insights into its variability throughout Earth's history. Community Climate System Model (CCSM4) simulations published in Ladant et al. (2020). CRET6x was then run for an additional ~ 1200 model years. The fully-coupled 3x PI CO2 Eocene simulation (EECO3x) has paleogeography, land-sea mask, and vegetation distribution following the Deep-Time Modeling Intercomparison Project (DeepMIP) protocol at ~ 55 million years ago (Herold et al. 2014). The simulation was run for ~ 2000 model years, and the ocean temperature and salinity were initialized from a Paleocene-Eocene Thermal Maximum (PETM) quasi-equilibrated state (Zhu et al. 2020). The fully-coupled 2x PI CO2 Oligocene simulation (OLIG2x) was run using iCESM1.2 for ~ 1600 model years. The paleogeography and land-sea mask used has been published by Straume et al. (2020). The vegetation distribution used is based on a global BIOME4 model for the Chattian (26 million years ago) by the Bristol Research Initiative for the Dynamic Global Environment (BRIDGE 2022). 2.3 Slab Ocean Model simulations We created SOM simulations using the climatological data from each equilibrated paleoclimate simulated ocean. A slab ocean is a simplified representation of the ocean, approximating a well-mixed ocean mixed layer at each grid point; SOM configurations have been traditionally used particularly for understanding the climate sensitivity of an environment through quicker, cost-effective means (Danabasoglu and Gent 2009; Singh et al. 2022; Zhu et al. 2019). SOM simulations use climatological data from fully-coupled, equilibrated oceans to prescribe spatially and seasonally varying ocean heat flux convergence (Qflux). The model includes a thermodynamic slab ocean with a prescribed mixed layer depth that varies spatially but does not simulate ocean currents. Surface salinity and temperature evolve prognostically, but there is no explicit advection or ocean dynamics (Kiehl et al. 2006). The Qflux is calculated using the density of seawater, ocean heat capacity, mixed layer depth, net heat flux into the ocean, and sea surface temperature (SST) (Bitz et al. 2012). A SOM simulation is still coupled with the other CESM components, but it only takes a few decades to equilibrate, even with an exceptionally high atmospheric CO2 level, and represents the foremost climatic processes that affect ECS (Danabasoglu and Gent 2009). We ran five SOM simulations, each using one of the distinct slab oceans and geographic configurations (CRET3x, CRET6x, EECO3x, OLIG2x, and PI1x) for ~ 100 years each and produced climatological files at that original CO2 level. We then doubled the CO2 level of each simulation and ran them for another ~ 100 years to produce new climatological files at the doubled atmospheric CO2 level. For the CRET3x and PI1x simulations, we doubled CO2 a subsequent time to achieve a second ECS at a higher CO2 background state. The OHT experiment was conducted by cloning the CRET6x SOM simulation, swapping its ocean climatological data with the cooler CRET3x ocean, running the simulation for ~ 100 years, doubling CO2, and running it again (see “ECS sensitivity to ocean heat transport”). The atmospheric and land conditions of each simulation in the OHT experiment are identical to those of the fully-coupled paleoclimate simulation from which the ocean climatology was derived. Page 5/23 = 3x PI CO2), and the last number represents the atmospheric CO2. The final simulation listed is for the OHT experiment (see “ECS sensitivity to ocean heat transport”). The rightmost columns describe whattotal percentage of the planet was covered by ocean in the initial SOM simulation, as well as specificallyhigh latitude ocean. Simulation Name ECS InitialGMST Initial %Ocean Initial % High Latitude (> 60°) Ocean CRET3xSOM3x◊CRET3xSOM6x 5.57 25.57 ℃ 73.0% 40.5% CRET3xSOM6x◊CRET3xSOM12x6.14 31.14 ℃ 73.1% 41.3% CRET6xSOM6x◊CRET6xSOM12x6.66 30.71 ℃ 73.1% 41.4% EECO3xSOM3x◊EECO3xSOM6x6.57 25.24 ℃ 73.6% 48.1% OLIG2xSOM2x◊OLIG2xSOM4x5.06 21.95 ℃ 72.1% 46.9% PI1xSOM1x◊PI1xSOM2x 4.04 14.73 ℃ 66.6% 22.0% PI1xSOM2x◊PI1xSOM4x 4.62 18.77 ℃ 68.6% 35.7% CRET3xOHT6x◊CRET3xOHT12x 5.77 31.08 ℃ 73.1% 41.3% While doubling CO2 increases GMST across all simulations, ECS varies by more than 2.5°C (Table 1). The highest ECS of 6.66°C occurs in the Cretaceous 6x PI CO2 case while the lowest, 4.04°C, is in the PI 1x PI CO2 case. The Cretaceous and Eocene ‘greenhouse’ climate simulations with higher initial GMST exhibit higher ECS than the Oligocene and PI ‘icehouse’ simulations. However, initial GMST alone does not determine ECS: for example, the CRET3xSOM6x simulation has the highest initial GMST at 31.14°C, yet only the third highest ECS at 6.14°C (Table 1, Fig. S7). In our simulations, a

Q-flux mixed layer model### Q-flux (mixed layer) model We can calculate the total freshwater mass and heat fluxes into the ice/ocean component from a spin up run with fixed SST. Given a climatology of ocean mixed layer depths, we can calculate the implied ocean heat convergence at the base of the mixed layer. This will incorporate the actual ocean heat transports, but also a residual component related to any errors in the surface heat or mass fluxes. This ocean heat convergence can be input into a thermodynamic mixed layer to give a Q-flux model. Depending on whether ice advection is turned on (see[above](#part4_2)), the horizontal transport of the sea ice is or is not incorporated into the ocean heat convergence. A full sea ice thermodynamic calculation is performed for this model. Generally this model takes 20 to 30 years to come into thermal equilibrium with any change in forcing, but this can be changed by specifying the maximum mixed layer depth (smaller implies a faster equilibration). There is an option to allow diffusion into the deep ocean by replacing OCNML with ODEEP in the run deck. This requires a 10 year spin up of the Q-flux model in order to properly set the deeper