Last Glacial Maximum (LGM) climate forcing and ocean dynamical feedback and their implications for estimating climate sensitivity

Equilibrium climate sensitivity (ECS) has been directly estimated using reconstructions of past climates that are different than today's. A challenge to this approach is that temperature proxies integrate over the timescales of the fast feedback processes (e.g., changes in water vapor, snow, and clouds) that are captured in ECS as well as the slower feedback processes (e.g., changes in ice sheets and ocean circulation) that are not. A way around this issue is to treat the slow feedbacks as climate forcings and independently account for their impact on global temperature. Here we conduct a suite of Last Glacial Maximum (LGM) simulations using the Community Earth System Model version 1.2 (CESM1.2) to quantify the forcing and efficacy of land ice sheets (LISs) and greenhouse gases (GHGs) in order to estimate ECS. Our forcing and efficacy quantification adopts the effective radiative forcing (ERF) and adjustment framework and provides a complete accounting for the radiative, topographic, and dynamical impacts of LIS on surface temperatures. ERF and efficacy of LGM LIS are −3.2 W m⁻² and 1.1, respectively. The larger-than-unity efficacy is caused by the temperature changes over land and the Northern Hemisphere subtropical oceans which are relatively larger than those in response to a doubling of atmospheric CO₂. The subtropical sea-surface temperature (SST) response is linked to LIS-induced wind changes and feedbacks in ocean–atmosphere coupling and clouds. ERF and efficacy of LGM GHG are −2.8 W m⁻² and 0.9, respectively. The lower efficacy is primarily attributed to a smaller cloud feedback at colder temperatures. Our simulations further demonstrate that the direct ECS calculation using the forcing, efficacy, and temperature response in CESM1.2 overestimates the true value in the model by approximately 25 % due to the neglect of slow ocean dynamical feedback. This is supported by the greater cooling (6.8 ^∘C) in a fully coupled LGM simulation than that (5.3 ^∘C) in a slab ocean model simulation with ocean dynamics disabled. The majority (67 %) of the ocean dynamical feedback is attributed to dynamical changes in the Southern Ocean, where interactions between upper-ocean stratification, heat transport, and sea-ice cover are found to amplify the LGM cooling. Our study demonstrates the value of climate models in the quantification of climate forcings and the ocean dynamical feedback, which is necessary for an accurate direct ECS estimation.

1 Introduction

Equilibrium climate sensitivity (ECS) is defined as the global mean surface air temperature (GMST) response to a doubling of atmospheric CO₂ and accounts for the Planck response and water vapor, ice albedo, lapse rate, and cloud feedbacks (with timescales < 100 years; “Charney Sensitivity”; Charney et al., 1979). By tradition and for practical reasons, ECS does not account for slow feedback processes, such as changes in vegetation, cryosphere and ocean circulation, effects of which are included in Earth system sensitivity (e.g., Lunt et al., 2010). The ECS range is estimated to be 2.6–3.9 ^∘C (66 % confidence interval) in a recent assessment (Sherwood et al., 2020), which represents a narrower range than the traditional one of 1.5–4.5 ^∘C (IPCC, 2013). Nevertheless, ECS is difficulty to constrain from present-day observations due to the brevity of the instrumental record, the small magnitude of radiative forcing and temperature response relative to the natural variability, and the dependence of ECS estimates on the transient sea-surface temperature (SST) pattern in the historical period (Knutti et al., 2017; Sherwood et al., 2020).

Paleoclimate records overcome these limitations and provide unique observational constraints on ECS. The Last Glacial Maximum (LGM; ∼21 ka) has been considered to be an ideal target for estimating ECS since it represents a quasi-equilibrium climate state with large changes in climate forcing and response and relatively high spatial coverage of well-dated proxy temperatures. PALAEOSENS Project Members (2012) proposed a framework to obtain ECS from reconstructions of paleo-temperatures and climatic forcings in which slow (timescales > 100 years) feedback processes such as changes in greenhouse gases (GHGs), land ice sheets (LISs), Earth's orbits, and land use are regarded as climate forcings rather than feedbacks. This approach has been widely used to directly calculate ECS from proxy reconstructions of the LGM and glacial–interglacial cycles with estimates ranging from 2.6 to 8.1 ^∘C (Friedrich et al., 2016; Köhler et al., 2017; Stap et al., 2019; Tierney et al., 2020; von der Heydt et al., 2014).

Direct ECS estimation using the PALAEOSENS Project Members (2012) approach relies on a complete understanding of the slow feedback processes. In the context of the LGM, both the radiative forcing due to changes in GHGs, LISs, vegetation, and aerosols and their efficacy (the ratio of the warming effect attributed to a given forcing relative to that due to a doubling of atmospheric CO₂ under preindustrial conditions) must be known. Of these forcings, GHGs and LISs have been considered in most LGM-based ECS estimations (Friedrich et al., 2016; Köhler et al., 2017; Schmittner et al., 2011; Stap et al., 2019; Tierney et al., 2020; von der Heydt et al., 2014). Previous estimates of the LGM LIS forcing account for albedo changes associated with the presence of LISs and the exposure of shelves due to the lowered sea level (Fig. 1a), yielding a shortwave forcing ranging from −1.5 to −5.2 W m⁻² (Braconnot et al., 2012; Braconnot and Kageyama, 2015; Friedrich et al., 2016; Hansen et al., 2013; Köhler et al., 2010; Taylor et al., 2007; Tierney et al., 2020). However, these estimates neglect changes in surface topography (Fig. 1b), which can change surface temperature and longwave emission. Moreover, LIS topographic changes altered atmospheric (Fig. 1c) and ocean circulations (Herrington and Poulsen, 2012; Kutzbach and Guetter, 1986; Zhu et al., 2014), which can change surface temperatures without directly involving radiative processes (e.g., through a wind–evaporation–SST feedback; Xie and Philander, 1994). To our knowledge, a complete quantification of the LGM LIS forcing that accounts for the radiative, topographic, and dynamical effects has not been done.

Figure 1 (a) Changes in shortwave surface albedo associated with the presence of the LGM ice sheets and shelf exposures due to the lowered sea level. Albedo changes are diagnosed using “fixed-SST” experiments (ATM_LIS and ATM_PI; see Sect. 2.2). Note the uneven color bar. (b) Changes in surface elevation associated with the LGM ice sheets. (c) Wind changes at 850 hPa as an illustration of the ice-sheet dynamical forcing. Shown are anomalies in the fixed-SST experiments (ATM_LIS and ATM_PI).

Furthermore, the efficacy of LIS forcing has received relatively little attention (Yoshimori et al., 2009), although it is clear that albedo effects are mostly distributed over high-latitude land. This oversight is problematic since the assumptions of LIS efficacy can greatly impact the resulting ECS estimates (Stap et al., 2019; Tierney et al., 2020).

Another caveat with the direct ECS estimation approach is that fast feedbacks depend on the climate state and interact with slow feedbacks. Global climate model (GCM) studies have shown that the magnitude of fast climate feedbacks varies with global temperature, e.g., stronger cloud and water vapor feedbacks at higher GMSTs (Caballero and Huber, 2013; Schneider et al., 2019; Yoshimori et al., 2009; Zhu and Poulsen, 2020b; Zhu et al., 2019). Furthermore, the LGM ocean circulation was different than today (e.g., Curry and Oppo, 2005). Ocean dynamical processes can influence global temperature through interactions with sea ice, SST pattern, and cloud processes (Dong et al., 2019; Ferrari et al., 2014; Rose et al., 2014; Shin et al., 2003; Winton et al., 2013; Zhou et al., 2017), constituting an ocean dynamical feedback that takes place on timescales longer than 100 years. The contribution of the ocean dynamical feedback to the magnitude of LGM cooling and its impact on direct ECS estimation have not been thoroughly studied.

In this study, we address whether ECS can be accurately estimated using the direct calculation approach and knowledge of the LGM climate forcing and global temperature. To answer this question, we adopt the adjusted forcing–feedback framework (Sherwood et al., 2015) to provide a complete quantification of the forcing and efficacy of LGM LIS and GHG using a suite of climate simulations, in comparison to previous studies that only considered surface albedo effects of LIS. We also investigate the role of the ocean dynamical feedback in modulating the magnitude and spatial distribution of the LGM cooling by comparing fully coupled and slab ocean simulations. Finally, we discuss the implications of our results for direct ECS estimation using paleoclimate reconstructions.

2 Method, model, and experiments

2.1 Model and fully coupled simulations

We employ the Community Earth System Model (CESM) version 1.2 with a horizontal resolution of 1.9×2.5^∘ (latitude × longitude) for the atmosphere and land and a nominal 1^∘ for the sea ice and ocean (Hurrell et al., 2013). CESM1.2 is among the models that best reproduce climate features from instrumental records (Knutti et al., 2013) and has been extensively used for studying past climates (DiNezio et al., 2016; Otto-Bliesner et al., 2015; Zhu et al., 2017, 2019, 2020). Our CESM1.2 experiments were run with prescribed satellite phenology (SP) in the land model (Community Land Model version 4; CLM4) without an active carbon–nitrogen (CN) biogeochemical cycle. In the SP mode, leaf area and stem area indices and vegetation heights in CLM4 are prescribed according to data derived from satellite observations (Lawrence et al., 2011). Our choice of CLM4 with satellite phenology is based on the overall poorer simulation of vegetation phenology with an active CN, which could potentially be more problematic for paleoclimate simulations (Lawrence et al., 2011).

We make use of the fully coupled preindustrial (PI) and LGM simulations (FCM_PI and FCM_LGM; Table 1) from our previous studies (Tierney et al., 2020; Zhu et al., 2017). FCM_LGM was forced by boundary conditions consistent with protocols from the Paleoclimate Modelling Intercomparison Project phase 4 (PMIP4), including altered GHG concentrations, Earth orbital parameters, and LISs (Kageyama et al., 2017). LIS forcing is derived from the ICE-6G reconstruction (Peltier et al., 2015) and includes changes in land elevation and surface properties due to the presence of LGM ice sheets as well as changes in the land–sea mask to account for the lower sea level. FCM_LGM used prescribed preindustrial vegetation cover and aerosol emissions, as reliable global reconstructions are not available (Kageyama et al., 2017; Köhler et al., 2010). Although FCM_LGM contains the orbital forcing, its effect on GMST is small (e.g., Liu et al., 2014) and neglected in the following analysis. The FCM_LGM ocean state was initialized from the LGM simulation using the Community Climate System Model version 4, which had been spun-up for more than 2400 years and reached a quasi-equilibrated state (Brady et al., 2013). FCM_LGM was integrated for an additional ∼1800 years to reach equilibration under an updated atmosphere model (DiNezio et al., 2016; Tierney et al., 2020; Zhu et al., 2017). The top-of-atmosphere (TOA) energy imbalance averaged over the last 100 years is −0.06 W m⁻² in FCM_LGM, which is comparable to the 0.09 W m⁻² in FCM_PI, indicating the surface climate has reached a quasi-equilibrium glacial state. The global volume-mean ocean temperature exhibited a cooling of 0.15 ^∘C in the last 900 years of FCM_LGM. GMST in FCM_LGM is 6.8 ^∘C lower than that in FCM_PI and falls within the range directly estimated from proxy data in a recent study (−6.8 to −4.4 ^∘C; Tierney et al., 2020).

Table 1 List of CESM simulations conducted in this study, including experiment name, climate forcing configurations, run length (years), GMST and GMST changes (^∘C), effective radiative forcing (ERF; W m⁻²), and efficacy (unitless). One standard deviation calculated using annual data is listed. Note that kernel-corrected ERF uses 12-month climatology and no uncertainty is provided.

2.2 “Fixed-SST” simulations and the effective radiative forcing

We adopt the forcing–feedback framework with the concept of rapid adjustments (Sherwood et al., 2015). We use fixed-SST experiments to calculate the effective radiative forcing (ERF), defined as the change in net TOA radiative flux after adjustments of the atmospheric temperature profile, water vapor, and clouds (Hansen et al., 2005). Our results show that this method is especially well-suited for quantifying the LIS forcing and is an advancement over either simplified bulk calculations or the approximate partial radiative perturbation method used in previous studies, which only provide an estimation of the shortwave forcing from albedo effects. In the fixed-SST experiments, an LGM climate forcing (e.g., GHG) is introduced into a preindustrial simulation with active atmosphere and land models, but with SST and sea ice prescribed to the unperturbed preindustrial climatology. The land surface temperature is allowed to adjust as it is impractical to fix in the model. The ERF attributed to a forcing is obtained as the change in TOA net radiation between simulations with and without the forcing. ERF and land temperature change in the fixed-SST experiments are termed ERF_fsst and ΔT_fsst, respectively. Changes in atmospheric temperature, water vapor, and clouds in response to the climate forcing, without mediation by the global-mean temperature, are referred to as “adjustments” (see Fig. 1 of Sherwood et al. (2015) for an illustration).

To quantify the ERF due to LGM LIS and GHG, we performed two fixed-SST experiments with LGM GHG and LIS forcing, respectively (ATM_GHG and ATM_LIS; Table 1). To examine whether ERFs from GHG and LIS are additive, we performed an additional experiment with both forcings included (ATM_LGM). To compare the forcing and response of LGM GHG and LIS to CO₂ increasing, we also carried out an experiment with twice the preindustrial atmospheric CO₂ concentration (ATM_2CO2). Finally, a standard preindustrial atmosphere-only simulation was performed and used as a reference (ATM_PI). These fixed-SST experiments used the same set of preindustrial SST and sea-ice coverage derived from the FCM_PI climatology. All fixed-SST simulations were run for 30 years with the last 25 used for analysis.

ERF_fsst contains radiation changes (biases) resulting from land surface temperature changes (ΔT_fsst). For example, ΔT_fsst is negative in ATM_GHG (see Fig. 2a), leading to an underestimation of the magnitude of ERF due to the decrease in Planck emission at lower land surface temperatures. To account for the radiative effects from ΔT_fsst, two corrected versions of ERF were computed. In the first, the ΔT_fsst effect on TOA radiation was corrected using the climate sensitivity parameter (α; units: K W⁻¹ m⁻²) as

$$\begin{array}{}\text{(1)}& {\mathrm{ERF}}_{\mathit{\alpha }}={\mathrm{ERF}}_{\mathrm{fsst}}-\mathrm{\Delta }{T}_{\mathrm{fsst}}/\mathit{\alpha }.\end{array}$$

α was obtained from the coupled simulation in a slab ocean configuration (see Sect. 2.3). In the second correction, radiative kernels were used:

$$\begin{array}{}\text{(2)}& {\mathrm{ERF}}_{\mathrm{kernel}}={\mathrm{ERF}}_{\mathrm{fsst}}-A_{\mathrm{Ts}}-A_{\mathrm{Ta}}-A_q-A_{\mathrm{alb}}.\end{array}$$

In this approach, ERF_kernel is obtained by subtracting the direct rapid adjustments associated with ΔT_fsst from ERF_fsst while keeping the indirect rapid adjustments, such as cloud responses (Tang et al., 2019). The direct rapid adjustments that are subtracted include effects over land from changes in surface temperature (A_Ts), tropospheric air temperature (A_Ta), tropospheric water vapor (A_q), and albedo (A_alb). A_Ta is calculated by assuming a constant lapse rate in the troposphere, i.e., the same tropospheric air temperature change as the surface. Similarly, A_q accounts for the effect from tropospheric water vapor change under the assumption of a constant lapse rate. A_q is calculated by scaling the total water vapor effect with the ratio between the temperature-induced radiative flux change from a constant lapse rate and that from the full tropospheric temperature change.

Figure 2 (a) Land surface temperature changes (ΔT_fsst) in response to the LGM GHG forcing diagnosed in “fixed-SST” experiment ATM_GHG. (b) ΔT_fsst in response to the LGM LIS forcing diagnosed in ATM_LIS. (c) The shortwave component of effective radiative forcing associated with LGM GHG forcing diagnosed in ATM_GHG. (e) As (c) but for the longwave component. (d, f) As (c, e) but for the LGM LIS forcing. Note that a present-day land–sea mask is shown in all the figures, which differs from the LGM mask due to a lower sea level; this will result in a temperature difference over the shelf exposure regions in (b). Surface temperature above sea ice is allowed to adjust in fixed-SST experiments.

2.3 Slab ocean model simulations and the efficacy of forcing

To compare temperature responses to different climate forcings and to estimate α, we performed sensitivity experiments in a slab ocean model (SOM) configuration without ocean dynamics (Table 1). The SOM uses prescribed mixed-layer depth and heat transport convergence (“q flux” hereafter) (Bitz et al., 2011) that are derived from FCM_PI. As a result, temperature responses in SOM simulations are caused by the fast feedback processes and exclude the ocean dynamical feedback. SOM_GHG includes LGM GHG levels and the preindustrial values of all the other boundary conditions. Similarly, SOM_LIS incorporates the LGM ice-sheet forcing including a higher topography, an altered land–sea distribution to account for effects from sea-level change, and modified land surface properties over ice sheets. To examine whether the climate responses are additive, we performed SOM_LGM, in which both LGM GHG and LIS forcings were added. In addition, we conducted SOM_2CO2 with CO₂ level 2 times the preindustrial value. These SOM simulations were integrated for 60 years to allow the model to reach equilibrium (with TOA energy imbalance < 0.1 W m⁻²) (Danabasoglu and Gent, 2009). Averages over the last 20 years are used for analysis.

The climate sensitivity parameter (α) was obtained from fixed-SST and SOM simulations and used to calculate ERF_α (Eq. 1). Specifically, for a climate forcing (GHG, LIS, or 2×CO₂), α was estimated as

$$\begin{array}{}\text{(3)}& \mathit{\alpha }={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{SOM}}-\mathrm{\Delta }{T}_{\mathrm{fsst}}}{{\mathrm{ERF}}_{\mathrm{fsst}}}}.\end{array}$$

ΔT_SOM is the equilibrated GMST response in the SOM simulation. ΔT_SOM−ΔT_fsst represents the SST-mediated surface air temperature changes that are associated with ERF_fsst.

We define the efficacy (ε) of a climate forcing (GHG or LIS) as a ratio of its temperature response scaled by its ERF_α to that of 2×CO₂:

$$\begin{array}{}\text{(4)}& \mathit{\epsilon }={{\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{SOM}}/{\mathrm{ERF}}_{\mathit{\alpha }}}{\left(\mathrm{\Delta }{T}_{\mathrm{SOM}}/{\mathrm{ERF}}_{\mathit{\alpha }}\right)}}}_{\mathrm{2}{\mathrm{CO}}_{\mathrm{2}}}.\end{array}$$

We note that the LGM GHG forcing and efficacy in this study are calculated using a “low-top” atmosphere model with prescribed stratospheric chemical tracers and exclude indirect effects from stratosphere chemistry (Hansen et al., 2005).

2.4 The radiative kernels and approximate partial radiative perturbation approach

To correct the ERFs of doubling CO₂ and LGM GHG and to understand their efficacy, we employ the radiative kernels that are developed for CESM (Pendergrass et al., 2018). In the analysis, we calculate changes in the 12-month climatology of the variable of interest (e.g., surface temperature) and multiply that by the corresponding radiative kernel to estimate the TOA radiation changes. Climate feedback parameters are obtained by normalizing the TOA radiation anomalies by the GMST changes. Kernels analyses are not performed for LGM LIS simulations, as the present-day kernels are not suitable due to the large difference in the characteristics of the forcing and response (Yoshimori et al., 2011).

The approximate partial radiative perturbation (APRP; Taylor et al., 2007) is used to quantify the shortwave forcing and feedback, in particular for the LIS simulations. In contrast to the radiative-kernel method, APRP is independent of the forcing and background climate state and produces results that differ from partial radiative perturbation (PRP) by less than 7 % (Taylor et al., 2007; Yoshimori et al., 2011; Zhu and Poulsen, 2020b; Zhu et al., 2019). APRP represents the atmosphere as a single layer with bulk optical properties and usually uses monthly mean model output to derive the radiative effects and feedbacks associated with changes in surface albedo, clear-sky processes, and clouds. The shortwave cloud feedback is further decomposed into contributions from changes in cloud amount, scattering, and absorption. APRP has been used in many previous studies to quantify the shortwave forcing associated with LIS albedo changes (Braconnot et al., 2012; Braconnot and Kageyama, 2015; Brady et al., 2013; Taylor et al., 2007; Tierney et al., 2020). When using APRP to quantify shortwave feedbacks, we only show results over model grid points that are ocean in both PI and LGM simulations; in this way, climate feedbacks are separated from forcings, e.g., from LIS, although feedback processes over land are overlooked.

3 Results

3.1 Effective radiative forcing

The ERF_fsst due to LGM GHG is $-\mathrm{2.6}\pm \mathrm{0.2}$ W m⁻² (±1σ; Table 1; Fig. 2c and e). After correcting the radiative effects associated with land temperature changes ($\mathrm{\Delta }{T}_{\mathrm{fsst}}=-\mathrm{0.2}$ ^∘C; Fig. 2a) using the climate sensitivity parameter and radiative kernels, the ERF_α  and ERF_kernel are $-\mathrm{2.8}\pm \mathrm{0.3}$ and −2.8 W m⁻², respectively, and agree well with previous estimates of −2.8 to −3.0 W m⁻² (Hansen et al., 2013; Köhler et al., 2010). For a doubling of CO₂, ERF_fsst, ERF_α, and ERF_kernel are 3.7±0.3, 3.9±0.3, and 4.0 W m⁻², respectively, well within the multi-model range in recent studies (Smith et al., 2018; Tang et al., 2019). For both LGM GHG and 2×CO₂, ERF_kernel falls in the middle of the uncertainty range of ERF_α, suggesting that both the correction methods using radiative kernels and climate sensitivity parameters produce meaningful and accurate results.

In response to LGM LIS, ΔT_fsst has a global mean of −1.3 ^∘C with maximum cooling over ice sheets exceeding −24 ^∘C, much greater than the land temperature changes associated with GHG forcing (Table 1; Fig. 2a and b). ΔT_fsst results from the higher surface albedo over regions with increased coverage of ice sheets and land (due to shelf exposures), the elevated ice-sheet topography, and radiative and dynamic atmospheric adjustments. The global mean ERF_fsst due to LGM LIS is $-\mathrm{1.9}\pm \mathrm{0.2}$ W m⁻², resulting from a shortwave component of −3.7 W m⁻² (ERF_{fsst_sw}; Fig. 2d) and a longwave component of 1.8 W m⁻² (ERF_{fsst_lw}; Fig. 2f). ERF_{fsst_sw} is lowest over ice sheets with values of less than −80 W m⁻². Using shortwave APRP, we attribute 77 % (−2.8 W m⁻²) of ERF_{fsst_sw} to surface albedo changes over regions of ice sheets and shelf exposure, 13 % (−0.5 W m⁻²) to surface albedo changes associated with snow cover increases outside the ice-sheet regions, and 8 % (−0.3 W m⁻²) to cloud adjustments. The majority of the cloud adjustments (−0.2 W m⁻²) occurs over the Indo-Pacific warm pool, where the exposure of the Sunda and Sahul shelves produces a surface cooling and drying (DiNezio et al., 2016) and increases cloud condensates through enhanced large-scale moist advection (Zhang et al., 2003). Outside the tropics, clouds diminish over ice sheets and in the downwind regions and shift with the position of the storm tracks; yet, the overall impact on the global mean ERF_{fsst_sw} is small. ERF_{fsst_lw} exceeds 40 W m⁻² over ice sheets, which results primarily from the reduced longwave radiation due to a higher effective emission elevation and lower temperatures (Figs. 1b and 2b).

Using the climate sensitivity parameter (diagnosed in SOM_LIS; Table 1; Eqs. 1 and 3), we calculate a global mean ERF_α from LGM LIS of $-\mathrm{3.2}\pm \mathrm{0.2}$ W m⁻². Our calculation accounts for the radiative, topographic, and dynamic adjustments associated with LIS, in contrast to only the albedo effects considered in previous studies (Braconnot et al., 2012; Braconnot and Kageyama, 2015; Taylor et al., 2007; Tierney et al., 2020). Using the APRP approach as in these previous studies, we calculate a shortwave forcing of −2.8 W m⁻² from surface albedo changes over regions with new ice-sheet coverage and shelf exposures in ATM_LIS. We note that this APRP approach overestimates the shortwave radiative forcing that is attributable exclusively to changes in LIS extent, as it includes the radiative effect of snow increases over ice sheets (or regions with shelf exposure); the albedo of fresh snow is considerably larger (0.8–0.9 versus ∼0.6) (Cuffey and Paterson, 2010). The LIS-induced cooling increases the proportion of snow relative to rain, which reflects more shortwave radiation than that from ice-sheet albedo alone. The snow-induced overestimation of the LIS forcing is larger if the cooling over ice sheets is greater. For example, the shortwave forcing from APRP analysis is greater in coupled simulations than that in atmosphere-only simulations with fixed PI SST (e.g., −3.3 W m⁻² in FCM_LGM versus −2.8 W m⁻² in ATM_LIS), due to the greater cooling over ice sheets.

Figure 3 Changes in surface temperature in (a) SOM_2CO2, (d) SOM_GHG, and (g) SOM_LIS. Changes in 850 hPa winds (units: m s⁻¹) in fixed-SST experiments are shown as vectors in (a) ATM_2CO2, (d) ATM_GHG, and (g) ATM_LIS. Changes in SST scaled by the corresponding GMST change in (b) SOM_2CO2, (e) SOM_GHG, and (h) SOM_LIS. The shortwave cloud feedback parameter over LGM ocean grid points diagnosed using the APRP approach in (c) SOM_2CO2, (f) SOM_GHG, and (i) SOM_LIS. (j) Zonal mean shortwave cloud feedback over ocean in SOM_2CO2 (red), SOM_GHG (blue), and SOM_LIS (green). Note that temperature and wind changes in 2×CO₂ experiments in (a) have been multiplied by −1 to facilitate the comparison with those in SOM_GHG and SOM_LIS.

3.2 Efficacy of LGM GHG and LIS forcings

Our results suggest that the efficiency of lowering GHGs to LGM levels is smaller than that of doubling atmospheric CO₂ under PI conditions, i.e., the LGM GHG forcing has a smaller-than-unity efficacy of 0.9±0.1 (Table 1; Eq. 4). In SOM_2CO2, GMST increases by 3.6 K in response to an ERF_α of 3.9 W m⁻². In SOM_GHG, GMST decreases by 2.2 K in response to an ERF_α of −2.8 W m⁻². The lower ε of LGM GHG forcing is caused by a weaker cloud feedback in response to cooling (Table 2). Using radiative kernels, we find that the Planck, albedo, and combined lapse rate and water vapor feedbacks stay largely unchanged; however, the cloud feedback parameter is 30 % smaller in SOM_GHG than in SOM_2CO2 (0.32 versus 0.46 W m⁻² K⁻¹). The decrease in the cloud feedback is due to the shortwave component; the cloud scattering feedback is weaker in response to cooling than that to warming over high-latitude regions, leading to a weaker shortwave response (Fig. 3c, f and j; Table 2 APRP columns). These results demonstrate a state-dependent cloud feedback that increases with GMST, a feature that has been found in the latest three CAM versions (Zhu and Poulsen, 2020b; Zhu et al., 2019) and many other climate modes (e.g., Crucifix, 2006).

Table 2 Climate feedback parameters (units: W m⁻² K⁻¹) in the SOM simulations. Radiative-kernel-based analysis is performed for SOM_2CO2 and SOM_GHG but not SOM_LIS due to the drastically different boundary conditions. APRP-based analysis is performed for all three simulations. APRP quantifies the shortwave climate feedback parameter and decomposes the cloud feedback into contributions from changes in cloud amount, scattering, and absorption. The cloud absorption feedback is not shown, as it is small and varies little between simulations. Values in parentheses are the contribution to the global mean value from ocean grid points in LGM simulation. Note that the high value in the albedo column for SOM_LIS (1.31) includes the contribution from the shortwave forcing over land.

The ε of LGM LIS is 1.1±0.1, resulting from an ERF_α of −3.2 W m⁻² and a ΔGMST of −3.2 K in SOM_LIS (Table 1) and suggesting that LIS forcing is 10 % more effective in changing GMST than doubling CO₂ when only fast feedbacks are considered. Overall, 40 % of the LIS-induced cooling (1.3 of 3.2 K) is attributed to land temperature changes that involve radiative, topographic, and dynamical effects of the LIS forcing and are independent from SST changes (Figs. 1 and 2; Table 1); land only accounts for 7 %–8 % of the GMST change in SOM_2CO2 and SOM_GHG. In addition to the large contribution from processes over land, the shortwave cloud feedback over ocean is greater in response to the LIS forcing (0.30 W m⁻² K⁻¹) than that to the doubling CO₂ forcing (0.21 W m⁻² K⁻¹; APRP analysis in Table 2 and Fig. 3). The greater shortwave cloud feedback is due to changes in both cloud amount and scattering and is especially prominent over the Northern Hemisphere subtropics (Fig. 3j), which likely reflects the southward shift of storm tracks and clouds. The process can also be understood as an “SST pattern effect”. The cloud feedback parameter is expressed as

$$\begin{array}{}\text{(5)}& {\mathit{\lambda }}_{\mathrm{CLD}}\equiv {\displaystyle \frac{\mathrm{dCRE}}{\mathrm{dGMST}}}={\displaystyle \frac{\partial \mathrm{CRE}}{\partial {\mathrm{SST}}_{\mathrm{SUB}}}}{\displaystyle \frac{{\mathrm{dSST}}_{\mathrm{SUB}}}{\mathrm{dGMST}}},\end{array}$$

where SST_SUB is the SST over the subtropical North Pacific and North Atlantic, and CRE denotes the global cloud radiative effects. SST_SUB is positively correlated with global CRE ($\frac{\partial \mathrm{CRE}}{\partial {\mathrm{SST}}_{\mathrm{SUB}}}>\mathrm{0}$) (Dong et al., 2019; Zhou et al., 2017). As a result, a greater SST_SUB change relative to the GMST gives rise to greater cloud feedback through changing the lower-tropospheric stability (Wood and Bretherton, 2006). Over the subtropical North Pacific and North Atlantic, the shortwave cloud feedback exceeds 10 W m⁻² K⁻¹ in SOM_LIS, in comparison to a maximum of ∼3 W m⁻² K⁻¹ in SOM_2CO2 (Fig. 3c and f), which is consistent with the relative magnitude of SST change in each experiment (Fig. 3b and h).

The formation of the SST responses over the subtropical North Pacific and North Atlantic is attributed to the ice-sheet-driven wind changes (Fig. 3g). In response to the topographic effects of LGM LIS, the Northern Hemisphere westerly jet shifts southward in ATM_LIS, producing cyclonic low-level wind anomalies over the subtropical and midlatitude North Pacific and anticyclonic anomalies over the subtropical North Atlantic (Kutzbach and Guetter, 1986; Zhu et al., 2014). This anomalous wind pattern forces regional SST changes through changing latent heat flux and amplify the coupled response through the wind–evaporation–SST feedback (Chiang and Bitz, 2005; Xie and Philander, 1994). For example, the trade wind strengthens over the subtropical North Atlantic, which cools the subtropical SST due to the enhanced evaporation and reinforces the anomalous wind pattern, forming a positive feedback. This non-radiative pathway of LIS's influence on the surface temperature is largely absent when GHGs are changed (Fig. 3a and d), highlighting the complex nature of non-GHG climate forcings and the importance of using efficacy to evaluate the overall effectiveness of their radiative forcing as compared to a doubling of atmospheric CO₂.

3.3 Are forcings and responses additive?

Our simulations suggest that ERFs and surface temperature responses of LGM GHG and LIS are globally additive. In fixed-SST experiments with both the LGM GHG and LIS forcings (ATM_LGM), ERF_fsst is $-\mathrm{4.4}\pm \mathrm{0.3}$ W m⁻², approximately the sum of those in ATM_GHG and ATM_LIS ($-\mathrm{2.6}\pm \mathrm{0.2}$ and $-\mathrm{1.9}\pm \mathrm{0.2}$, respectively). Similarly, the ERF_α in ATM_LGM, $-\mathrm{6.1}\pm \mathrm{0.3}$ W m⁻², is nearly equal to the sum of those in ATM_GHG and ATM_LIS ($-\mathrm{2.8}\pm \mathrm{0.3}$ and $-\mathrm{3.2}\pm \mathrm{0.2}$ W m⁻², respectively). The SOM_LGM ΔGMST in response to combined GHG and LIS forcings is $-\mathrm{5.3}\pm \mathrm{0.09}$ ^∘C and is close to the sum of ΔGMSTs in SOM_GHG and SOM_LIS ($-\mathrm{2.2}\pm \mathrm{0.11}$ and $-\mathrm{3.2}\pm \mathrm{0.09}$ ^∘C, respectively). From these results, we conclude that ERF and ΔGMST due to individual forcings are additive at the global level, which supports the approach to separate the LGM climate forcing and response into components associated with individual forcing agents. We note that at the regional level, especially over high latitudes, the ERF and ΔGMST from the sum of individual forcings and combined forcings do not match as well as at the global level (figure not shown), likely due to local feedbacks related to sea ice. It remains unclear whether the temperature responses to individual forcings are additive when slow feedback processes (e.g., ocean dynamics) are included.

Figure 4 (a) LGM cooling in the fully coupled model simulations (FCM_LGM–FCM_PI). (b) LGM cooling in the slab ocean model simulations (SOM_LGM–SOM_PI). Both SOM_LGM and SOM_PI use the same ocean mixed-layer depth and heat flux convergence (Q_flx) that are derived from the fully coupled preindustrial simulation (FCM_PI). (c) Difference in the simulated LGM cooling between SOM and FCM (a, b). (d) As (c) but for the SOM simulation with the prescribed Q_flx over the Southern Ocean (90–40^∘ S) replaced with that from FCM_LGM. (e) As (c) but for the SOM simulation with the prescribed Q_flx over 90^∘ S–30^∘ N replaced with that from FCM_LGM. (f) As (c) but for the SOM simulation with the prescribed Q_flx over the global ocean replaced with that from FCM_LGM. Note that a small intrinsic bias in surface temperature associated with SOM simulations (e.g., SOM_PI–FCM_PI) has been subtracted when comparing SOM and FCM simulations.

3.4 The ocean dynamical feedback

Our results demonstrate that the full extent of LGM cooling cannot be produced using an SOM configuration that accounts for fast feedback processes but excludes the slow ocean dynamical changes. ΔGMST is −6.8 ^∘C in FCM_LGM and −5.3 ^∘C in SOM_LGM (Table 1; Fig. 4a and b). Both simulations have reached an equilibrium state at the surface under the same climate forcings and only differ in the complexity of the ocean model. The LGM cooling is 1.5 ^∘C (28 %) greater with active ocean dynamics and interactions with the atmosphere and sea ice. The difference in LGM cooling between SOM and FCM primarily occurs in the Southern Ocean (SO), where SOM_LGM simulates a weaker LGM cooling by more than 10 ^∘C (Fig. 4a–c). In the eastern equatorial Pacific and eastern subtropical oceans (except for the subtropical Atlantic), SOM_LGM simulates a smaller LGM cooling by 1–2 ^∘C. In the North Atlantic, LGM cooling in SOM_LGM is greater by ∼1 ^∘C in the subtropics and more than 5 ^∘C along the sea-ice margin in subpolar regions. Over the Indo-Pacific Warm Pool and the western subtropics, surface temperature change is similar between FCM and SOM simulations, suggesting a limited role of ocean dynamical response over these regions.

Accounting for the SO dynamical effects increases the LGM cooling in the SOM configuration and explains the majority (67 %) of the difference from FCM simulations. This is shown in an LGM SOM simulation (SOM_SO), in which the prescribed q flux over the SO (90–40^∘ S) is replaced with those derived from FCM_LGM, with other regions remaining unchanged (using values from FCM_PI). SOM_SO simulates a colder LGM than SOM_LGM, especially over the SO, where the large temperature difference (>10 ^∘C) between SOM_LGM and FCM_LGM is mostly removed (Fig. 4d versus Fig. 4c). In addition to the impact on local temperatures, SOM_SO simulates lower surface temperatures over the eastern equatorial Pacific and Indian Ocean and the Southern Hemisphere subtropics, producing a better match with FCM_LGM over these regions. This reflects a remote impact of the SO processes on the lower latitudes through changing tropical atmospheric circulations (Hwang et al., 2017).

The SO dynamical effects primarily result from upper-ocean stratification changes and the coupling with sea ice. In FCM_PI, the SO is stratified with the maximum ocean temperature occurring in the subsurface (500–1000 m; Fig. 5a). In zonal and annual mean, isotherms (shadings in Fig. 5a) intersect isopycnals (contours) mostly near 65–60^∘ S and along the Antarctic coast, indicating strong heat diffusion towards the mixed layer, i.e., a heat flux convergence of approximately −20 W m⁻² (Fig. 5c; red curve). The strong heat flux convergence warms the mixed layer and inhibits sea-ice formation, resulting in a quasi-permanent sea-ice extent to 68^∘ S (red horizontal bar in Fig. 5c; defined using a 70 % annual mean sea-ice cover). In comparison, the upper-ocean stratification in FCM_LGM is greatly reduced with a potential density change of less than 0.4 kg m⁻³ over most water columns and a largely invariant ocean temperature of −2 ^∘C (Fig. 5b). Meanwhile, the center of mixed-layer heat flux convergence is shifted northward to ∼56^∘ S and the quasi-permanent sea ice expands to ∼58^∘ S (Fig. 5c). An initial LGM cooling in the SO (e.g., caused by fast feedback processes) increases sea-ice formation and brine rejection, which enhances convection and decreases the upper-ocean stratification, resulting in a decrease in heat flux convergence to the mixed layer and amplifying sea-ice expansion. This feedback loop is absent in an SOM configuration with prescribed mixed-layer depth and heat flux convergence, leading to little expansion of the quasi-permanent sea ice (cyan horizontal bar in Fig. 5c) and much less LGM cooling in SOM_LGM (Fig. 4b and c). When the LGM changes in mixed-layer depth and heat flux convergence are prescribed in SOM_SO, sea ice expands northward (light green bar in Fig. 5c) and the SO experiences cooling (Fig. 4d versus Fig. 4c).

Figure 5 (a) Zonal mean potential temperature (shadings; units: ^∘C) and potential density (contours) over the Southern Ocean in the fully coupled preindustrial simulation (FCM_PI). Potential density values are reported in kilograms per cubic meter subtracting 1000 kg m⁻³. (b) As (a) but for the fully coupled LGM simulation (FCM_LGM). (c) Zonal mean heat flux convergence (units: W m⁻²) in the ocean mixed layer in FCM_PI (red) and FCM_LGM (blue). Negative values indicate that ocean dynamical processes transport heat from below into the mixed layer. Horizonal bars indicate the zonal mean semi-permanent sea-ice extent, defined as a 70 % of the annual mean sea-ice cover, in the FCM_PI (red), FCM_LGM (blue), SOM_LGM (cyan), and SOM simulations with Q_flx over 90–40^∘ S (light green) or 90^∘ S–30^∘ N (dark green) replaced with that in FCM_LGM.

Accounting for additional ocean dynamical effects in the low latitudes and the Northern Hemisphere further decreases the difference in LGM cooling between SOM and FCM simulations. This is supported by additional SOM simulations, in which we replace the prescribed q flux with those derived from FCM_LGM over 90^∘ S–30^∘ N and the entire global ocean. The tropical ocean dynamical effects decrease SST in the eastern equatorial Pacific and the Southern Hemisphere subtropics (Fig. 4e). Ocean dynamics in the Northern Hemisphere middle and high latitudes increases SST over the North Atlantic (Fig. 4f). The tropical ocean dynamical effects primarily result from changes in tropical ocean circulations and the coupling with the atmosphere (DiNezio et al., 2011; Vecchi and Soden, 2007). The Northern Hemisphere ocean dynamical effects are related to a stronger Atlantic Meridional Overturning Circulation (AMOC) and a greater northward ocean heat transport in the model (Brady et al., 2013). After accounting for the global ocean dynamical effects by using the q flux derived from FCM_LGM, high-latitude oceans still exhibit a temperature difference of ∼1–2 ^∘C, contributing to a GMST of ∼0.3 ^∘C, likely reflecting the challenge of using prescribed q flux to approximate the full extent of ocean dynamical effects (Bitz et al., 2011). Nevertheless, these results highlight the important role of the slow ocean dynamical feedback in modulating regional and global temperatures.

4 Discussion: implications for estimating climate sensitivity

Results presented herein highlight major caveats of the direct ECS estimation approach. Firstly, a complete understanding of the magnitude and efficacy of forcing agents is necessary, especially for non-GHG forcings (e.g., LIS, vegetation, and aerosols) that may have distinct spatial distribution and non-radiative pathways to change the energy balance of Earth. We suggest a GCM-based approach using the effective radiative forcing and adjustment framework to account for the complicated aspects of paleoclimate non-GHG forcings. In this approach, a fixed-SST simulation of ∼30 years with a forcing of interest is first conducted to calculate the effective radiative forcing (ERF_fsst) and the associated land temperature changes. An ERF_α is then obtained by correcting the ERF_fsst using the climate sensitivity parameter that is derived in an additional SOM simulation of ∼60 years. Moreover, efficacy of the forcing can also be derived using these simulations. This approach provides a complete consideration of the radiative and non-radiative effects of the forcing agent and is more consistent with the basic definition of the forcing–feedback framework. In contrast, the APRP-based approach used in previous studies only accounts for the effects from albedo changes. We note that, due to the inclusion of snow effects in the forcing quantification, the APRP-based approach overestimates the direct shortwave albedo effects that are attributable only to changing LIS extent. Our simulations suggest an LGM LIS efficacy of 1.1, which differs from the 0.45 in Stap et al. (2019). A precise explanation about this difference is challenging, given the large differences in the definition of forcing or efficacy, model complexity, and experimental design. Future studies are needed to quantify the effective radiative forcing and efficacy of the LGM vegetation and aerosols, as they are prescribed at the preindustrial values in our LGM simulations.

A second caveat concerns the role of the ocean dynamical feedback, which occurs on timescales of 10²–10³ years and should be accounted for when directly estimating ECS using forcing or response of an equilibrium climate. This complication stems from defining ECS to include only fast feedback processes with timescales of less than 100 years. Ocean feedback processes, including the heat redistribution by ocean circulation and the coupling with the atmosphere or sea ice, require more than 100 years to develop. Reconstructions of past climate forcings and GMST usually do not directly constrain ocean circulations and therefore could potentially impact the ECS estimation.

Figure 6 ECS in CESM1.2 (horizontal line) and the LGM-based estimates (box-and-whisker plots) using different values of forcing efficacy and ocean dynamical feedback (ODF)-induced ΔGMST changes. Each ECS estimation is obtained by performing 10 000 Monte Carlo calculations, which incorporates the uncertainties (assumed to be Normal) in forcings and temperature responses. The box and whisker indicate a 68 % and 95 % confidence interval, respectively. ERFs and efficacy of LGM GHG and LIS are listed in Table 1. ΔGMST changes from the ocean dynamical feedback is the difference between FCM and SOM simulations (Table 1).

To demonstrate the above caveats, we assume CESM1.2 is a perfect model and estimate ECS using LGM constraints that are derived from model simulations as

$$\begin{array}{}\text{(6)}& \mathrm{ECS}={\displaystyle \frac{\mathrm{\Delta }\mathrm{GMST}-\mathrm{\Delta }{\mathrm{GMST}}_{\mathrm{ODF}}}{{\mathit{\epsilon }}_{\mathrm{GHG}}\times {\mathrm{ERF}}_{\mathrm{GHG}}+{\mathit{\epsilon }}_{\mathrm{LIS}}\times {\mathrm{ERF}}_{\mathrm{LIS}}}}{\mathrm{ERF}}_{\mathrm{2}{\mathrm{CO}}_{\mathrm{2}}}.\end{array}$$

In the above equation, ΔGMST_ODF denotes the GMST change (approximately −1.5 ^∘C; see Sect. 3.4) that is caused by the slow ocean dynamical feedback and is subtracted from the total LGM cooling (ΔGMST = −6.8 ^∘C in CESM1.2). ERF_GHG, ERF_LIS, and ERF${}_{\mathrm{2}{\mathrm{CO}}_{\mathrm{2}}}$ in our simulations are $-\mathrm{2.8}\pm \mathrm{0.3}$, $-\mathrm{3.2}\pm \mathrm{0.2}$, and 3.9±0.3 W m⁻², respectively (Table 1). ε_GHG and ε_LIS are 0.9 and 1.1, respectively. In our “perfect model” assumption, all the above values are unbiased, and the “true” ECS is 3.6 ^∘C. We perform ECS calculations using Eq. (6) with 10 000 Monte Carlo draws to sample the uncertainty in forcings and explore impacts from different assumptions of climate forcing or efficacy and the ocean dynamical feedback (Fig. 6). If we do not remove the effects of the ocean dynamical feedback and assume that both GHG and LIS forcings have a unit efficacy, as has been done in most previous studies (e.g., PALAEOSENS Project Members, 2012), we obtain a median ECS of 4.5 ^∘C, an overestimate of 25 % that is statistically distinguishable from uncertainties associated with climate forcings. Using the true efficacy of the LGM GHG or LIS produces a small change in ECS ($<\sim \mathrm{0.3}$ ^∘C) that approximately cancels each other, as ε_GHG is smaller and ε_LIS is larger than unity. Accounting for the ocean dynamical feedback greatly improves the ECS calculation, yielding a median of 3.5 ^∘C, 0.1 ^∘C smaller than the true ECS. We note that here we account for the ocean dynamical effect by subtracting the corresponding contribution from the total LGM GMST change. Alternatively, we can also invoke a non-constant sea-ice albedo feedback that depend on ocean dynamics (see Fig. 5c). Nevertheless, this exercise highlights the importance of the ocean dynamical processes, which, if not accounted for, may cause an overestimation of the (“fast feedback”) ECS value using reconstructions of LGM forcings or responses.

The tight coupling between sea-ice extent and ocean dynamics in the Southern Ocean identified in our simulations is consistent with previous modeling and theoretical studies (Ferrari et al., 2014; Shin et al., 2003). The quantitative contribution of the ocean dynamical feedback to the LGM ΔGMST is likely model dependent, which, we speculate, could partly explain the lack of correlation between global and regional mean LGM cooling and ECS in Paleoclimate Modelling Intercomparison Project models (Hargreaves et al., 2012; Hopcroft and Valdes, 2015). Major features of the ocean circulation and seawater characteristics in our LGM simulations agree well with findings from proxy reconstructions (e.g., Adkins et al., 2002; Curry and Oppo, 2005), including an expansion of the Antarctic Bottom Water, a shallower North Atlantic Deep Water, an increase in abyssal stratification, and a saltier and colder southern-source deep water; however, a detailed examination of these features is beyond the scope of this study. Due to limited computing resources, our fully coupled LGM simulation has a cooling trend in the deep ocean (see Sect. 2.1), which will not impact our results on LGM radiative forcing or efficacy but will likely cause an underestimation of LGM ΔGMST and the importance of the ocean dynamical feedback in the model.

5 Conclusions

In this study, we have quantified the radiative forcing and efficacy of LGM GHGs and LISs in CESM1.2 and examined the contribution of the ocean dynamical feedback to surface temperature changes by comparing simulations in fully coupled and slab ocean configurations. ERFs of LGM GHG and LIS are estimated to be −2.8 and −3.2 W m⁻², respectively. The efficacy of LGM GHG and LIS forcings are estimated to be 0.9 and 1.1, respectively, indicating that lowering GHGs to LGM levels is 10 % less efficient in changing global temperature than that of doubling atmospheric CO₂ under PI conditions, while the LGM LIS is 10 % more efficient. The smaller-than-unity efficacy of LGM GHG forcing is primarily attributed to a smaller shortwave cloud feedback at lower temperatures, which is consistent with previous studies showing a temperature-dependent cloud feedback over high latitudes (e.g., Zhu and Poulsen, 2020b). The greater-than-unity efficacy of LGM LIS forcing is caused by relatively larger temperature changes over land and the Northern Hemisphere subtropical oceans, which are linked to the LIS-induced wind changes and feedbacks in ocean–atmosphere coupling and clouds. Our calculations of LIS forcing and efficacy account for the radiative effects from ice-sheet albedo and the topographic and dynamic effects associated with the ice-sheet elevation, in contrast to previous estimation that only considered the former. In addition, our simulations suggest that the effective radiative forcings and surface temperature responses of LGM GHG and LIS forcings are additive on the global level, which supports the approach in which individual forcing agents are considered separately.

Our simulations demonstrate that the full extent of LGM cooling cannot be realized if only fast feedbacks are accounted for. Overall, the slow ocean dynamical feedback amplifies the LGM cooling by 28 % (from 5.3 to 6.8 ^∘C) in CESM1.2. LGM-based ECS calculations that fail to account for this ocean dynamical effect produce an overestimation of fast feedbacks (by approximately 25 % in CESM1.2). In our simulations, the ocean dynamical feedback is primarily attributed to dynamical changes in the Southern Ocean, where a dynamical interaction between the upper-ocean stratification, mixed-layer heat flux convergence, and sea-ice cover is found to amplify the LGM cooling. Additionally, dynamical processes in the tropical oceans and the Atlantic also impact the regional and global temperatures. Overall, our results suggest an important role of climate models in the quantification of climate forcings and efficacy and the ocean dynamical feedback.