ZEMBA v1.0: an energy and moisture balance climate model to investigate Quaternary climate

The Zonally Averaged Energy and Moisture BAlance (ZEMBA) climate model is introduced as a simple and computationally efficient tool for studies of the glacial–interglacial cycles of the Quaternary. The model is based on an energy balance model comprising an atmospheric layer, a land component and a two-dimensional ocean transport model with sea ice. In addition, ZEMBA replaces temperature with moist static energy for calculations of diffusive heat transport in the atmospheric layer and includes a hydrological cycle for simulating precipitation and snowfall. Prior to coupling with an ice sheet model, we present and evaluate equilibrium simulations of the model for the pre-industrial period and the Last Glacial Maximum, using prescribed land ice fractions and elevation. In addition, we test the sensitivity of ZEMBA to a doubling of the atmospheric CO₂ concentration and a 2 % increase in solar radiation at the top of the atmosphere. Compared to a global climate model (the Norwegian Earth System Model version 2, NorESM2) and reanalysis data (ERA5), ZEMBA reproduces the zonally averaged climate of the pre-industrial period with reasonable accuracy, capturing features such as surface temperature, precipitation, radiative fluxes, snow cover, sea ice cover and meridional heat transport. The response of ZEMBA to increasing CO₂ concentrations is qualitatively similar to the observational record and climate models of higher complexity, including polar amplification over the Northern Hemisphere and during the winter months. The globally averaged rise in surface air temperature for a doubling in CO₂ is 3.6 °C. Finally, ZEMBA shows success in emulating changes in surface temperature and precipitation during the Last Glacial Maximum when compared to reconstructions and global climate models.

1 Introduction

Since the beginning of the Quaternary period (2.58 Ma to present), the Earth's climate has repeatedly switched between cold “glacial” periods and warmer “interglacial” periods, which are collectively known as glacial–interglacial cycles. Glacial periods are characterized by the presence of large ice sheets covering North America and Fennoscandia, while interglacials refer to times when ice sheets are restricted to Greenland and Antarctica, such as the current Holocene period (11.7 ka to present). The glacial–interglacial cycles are widely documented in natural climate archives. For example, the oxygen isotope ratio (δ¹⁸O) recorded in the shells of micro-organisms that accumulated on the ocean floor – hereafter referred to as benthic δ¹⁸O – serves as a valuable proxy for global ice volume and deep-ocean temperatures spanning millions of years into the past (Lisiecki and Raymo, 2005; Elderfield et al., 2012). Additionally, Antarctic ice cores contain valuable information over the past 800 kyr such as atmospheric CO₂ concentrations (Bereiter et al., 2015) and surface temperature changes (Jouzel et al., 2007; Kawamura et al., 2017).

Despite extensive research, a comprehensive understanding of what caused the glacial–interglacial cycles has remained elusive, although the importance of changes in the Earth's orbital parameters has been widely acknowledged. Changes in these orbital parameters are thought to determine when the Earth switches between glacial and interglacial climates by redistributing the incoming solar radiation (insolation) the Earth receives across latitudes and seasons. These variations in the orbital parameters encompass changes in the shape of the Earth's orbit (eccentricity) on 100 and 413 kyr cycles, changes in the tilt (obliquity) of the Earth's rotational axis on 41 kyr cycles, and changes in the time of the year when the Earth is closest to the Sun (climatic precession) on 19 and 23 kyr cycles. The most favoured of the orbital hypotheses comes from Milutin Milanković, who proposed that glaciation occurs during times of reduced insolation at the high northern latitudes during the summer months, when obliquity is low and the Northern Hemisphere (NH) summers coincide with the Earth's furthest distance from the Sun (and vice versa for deglaciation). Hays et al. (1976) were the first to discover that benthic δ¹⁸O records contained cycles with periods of 23, 41 and 100 kyr, in correspondence to the Earth's orbital cycles. Subsequent research has extended the benthic δ¹⁸O record further back in time (Pisias and Moore, 1981; Ruddiman et al., 1986) and led to compiled records from across the globe (Ahn et al., 2017; Lisiecki and Raymo, 2005), which together with other proxy records (Jouzel et al., 2007; Kawamura et al., 2017) provides compelling evidence for an orbital control on climate change. Nevertheless, despite a distinct orbital “rhythm” to glacial–interglacial cycles, important characteristics of the benthic δ¹⁸O record cannot be easily explained by Milanković theory, including a shift in the dominant periodicity observed around 1 Ma.

The Mid-Pleistocene Transition (MPT) from 1.25 to 0.7 Ma represents a transition in the dominant periodicity of glacial–interglacial cycles from 41 kyr in the Early Pleistocene to ∼ 100 kyr in the Late Pleistocene, occurring without significant variations in the orbital cycles (Clark et al., 2006). Hypotheses for the MPT are abundant in the literature (Berends et al., 2021). Prominent theories include the removal of a subglacial regolith beneath the NH ice sheets (Clark and Pollard, 1998; Clark et al., 2006), a gradual decline in atmospheric CO₂ concentrations (Raymo et al., 1988; Berger et al., 1999), the merging of the Laurentide and Cordilleran ice sheets over North America (Bintanja and Van De Wal, 2008), phase locking between Antarctica and the NH ice sheets (Raymo et al., 2006), or some combination of these mechanisms (Chalk et al., 2017; Willeit et al., 2019). Even the ∼ 100 kyr cycles in glacial ice volume following the MPT were not an expected outcome of Milanković theory, as they correspond to changes in the Earth's eccentricity, which has a negligible direct influence on summer insolation (Imbrie et al., 1993). Instead, these ∼ 100 kyr cycles have been explained as the skipping of one (80 kyr) or two (120 kyr) obliquity cycles (averaging to 100 kyr), leading to longer glacial cycles in the Late Pleistocene (Huybers and Wunsch, 2005), where the precise timing of deglaciation is set by the combined forcing of precession and obliquity (Huybers, 2011; Parrenin and Paillard, 2012). Independent modelling studies have also highlighted the importance of precession, and its modulation by eccentricity, for generating the ∼ 100 kyr cycles of the Late Pleistocene (Abe-Ouchi et al., 2013; Ganopolski and Calov, 2011). Interestingly, models that have simulated the MPT and/or the ∼ 100 kyr cycles of the Late Pleistocene still struggle to reproduce the dominant 41 kyr cycles seen in the Early Pleistocene (Berger et al., 1999; Willeit et al., 2019; Watanabe et al., 2023).

While modelling experiments of the Early Pleistocene capture 41 kyr cycles in ice volume, in correspondence to the Earth's obliquity, they generate stronger 19 and 23 kyr precession cycles (Berger et al., 1999; Willeit et al., 2019; Watanabe et al., 2023) than observed in the benthic δ¹⁸O record (Lisiecki and Raymo, 2005). This is unsurprising given that climatic precession controls the intensity of summer insolation. Indeed, for various metrics of summer insolation variability, including mid-month insolation (i.e. 21 June), monthly mean insolation or the caloric summer half year, precession has a strong influence at latitudes where NH ice sheets grow and melt. It is worth noting that precession cycles are detectable prior to the MPT (Lisiecki and Raymo, 2007; Liautaud et al., 2020), but they become more pronounced across the Quaternary, with much stronger signals observed in the Late Pleistocene compared to the Early Pleistocene (Raymo and Nisancioglu, 2003). Theories that account for weaker 19 and 23 kyr precession cycles during the Early Pleistocene include a counterbalancing between summer insolation intensity and summer duration (Huybers, 2006; Huybers and Tziperman, 2008), the cancellation of out-of-phase precession cycles between the NH ice sheets and Antarctica (Raymo et al., 2006; Morée et al., 2021), or the influence obliquity has on the poleward flux of moisture and accumulation rates on ice sheets (Raymo and Nisancioglu, 2003; Nisancioglu, 2004).

Models of varying levels of complexity have been employed to address these questions relating to glacial–interglacial cycles. Global climate models (GCMs) are too computationally expensive for simulations on these timescales, so studies have instead relied on conceptual models (Paillard, 1998; Parrenin and Paillard, 2003; Legrain et al., 2023), energy balance models (Pollard, 1978; Huybers and Tziperman, 2008), Earth system models of intermediate complexity (Ganopolski and Calov, 2011; Willeit et al., 2019), climate parameterization based on discrete GCM “snapshots” (Abe-Ouchi et al., 2013) or benthic δ¹⁸O records (Bintanja and Van De Wal, 2008). Among these, zonally averaged energy balance models (EBMs) provide computationally efficient tools for studying the response of the Earth's climate to changes in the Earth's orbital parameters (Suarez and Held, 1979; Pollard, 1978; Huybers and Tziperman, 2008; Stap et al., 2014). EBMs calculate the distribution of surface temperature with latitude by considering the conservation of energy on a sphere subject to heating by solar insolation, cooling by terrestrial radiation and the diffusive redistribution of heat from the Equator to the poles. Since the original works of Budyko (1969) and Sellers (1969), EBMs have long been used to study the Earth's climate sensitivity and have shown success at simulating both present-day and glacial climate states (North, 1975; Peng et al., 1987; Harvey, 1988; Jentsch, 1991; Bintanja and Oerlemans, 1996; Bintanja, 1997; Stap et al., 2014). Moreover, the simplicity of EBMs enables the isolation and identification of important processes and feedbacks (Huybers and Tziperman, 2008; Stap et al., 2014). Consequently, EBMs provide valuable tools for exploring the glacial–interglacial cycles of the Quaternary and guiding further investigations with more realistic models.

In this study, we introduce the Zonally Averaged Energy and Moisture BAlance (ZEMBA) climate model to study the response of the Earth's climate to changes in the orbital parameters. The model is designed to place a particular emphasis on physical processes that may influence the relative contributions of obliquity and precession to climate variability during periods such as the Early Pleistocene. The model spans both hemispheres and is forced by the full seasonal cycle in insolation. Unlike previous EBMs used for studies of glacial–interglacial cycles (Pollard, 1978; Huybers and Tziperman, 2008; Stap et al., 2014), ZEMBA includes a hydrological cycle to simulate precipitation and snowfall. Before using ZEMBA in experiments relating to glacial–interglacial cycles, it is important to ensure the model can simulate the present-day (or pre-industrial) climate with reasonable accuracy and to constrain the sensitivity of the model. Consequently, the following sections document the model and its climate sensitivity. Firstly, we describe a “control” simulation of ZEMBA for the pre-industrial period and compare the broad-scale features with a GCM and atmospheric reanalysis (Sect. 3.1). Then, we test the sensitivity of the model for a doubling of the atmospheric CO₂ concentration and a 2 % increase in solar insolation (Sect. 3.2). Finally, we evaluate its performance for a simulation of the Last Glacial Maximum (LGM) (Sect. 3.3). Prior to that, we provide a detailed description of the atmospheric (Sect. 2.1), land (Sect. 2.2) and ocean (Sect. 2.3) components of ZEMBA.

2 Model description

ZEMBA is primarily based on the EBM from Bintanja (1997), which comprises a single atmospheric layer overlying a surface divided into land and ocean. While utilizing the same shortwave and longwave radiation scheme as in Bintanja (1997), in addition to the same ocean transport model, ZEMBA includes a hydrological cycle to estimate precipitation and snowfall rates. Moreover, instead of parameterizing the surface albedo over land as a function of surface temperature, ZEMBA directly estimates snow coverage through the competition between snow accumulation (from the hydrological cycle) and ablation (from the surface energy balance). Following recent studies (Hwang and Frierson, 2010; Rose et al., 2014; Roe et al., 2015; Feldl and Merlis, 2021), atmospheric heat transport is now proportional to meridional gradients in near-surface moist static energy – instead of classic (or “dry”) EBMs that diffuse heat along temperature gradients. Moreover, a Hadley cell parameterization is included to produce an equatorward flux of moisture in the tropics (Siler et al., 2018). The radiative and turbulent heat fluxes are calculated separately over land and ocean, but atmospheric temperatures and humidities are set to the zonal average of land and ocean at the end of each model time step. The subsequent sections will describe in more detail the individual atmospheric, land and ocean components of ZEMBA.

2.1 Atmosphere

A vertically averaged atmospheric layer simulates the evolution of near-surface air temperature (T_a) through time. It evolves as a function of the radiative fluxes exchanged at the top and bottom of the atmospheric layer, the turbulent exchange of heat with the surface layer, the latent heat released during precipitation and snowfall, and the divergence of atmospheric heat transport. The temporal evolution of T_a is described as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{C}_{\mathrm{a}}{H}_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{a}\left(i\right)}}{\mathit{\delta }t}}=S_{\mathrm{a}\left(i\right)}+I_{\mathrm{a}\left(i\right)}+K_{\left(i\right)}+L_{\mathrm{v}}{P}_{\mathrm{r}\left(i\right)}+L_{\mathrm{f}}{P}_{\mathrm{s}\left(i\right)}\\ \text{(1)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}^{\mathrm{2}}\cos \left(\mathit{\theta }\right)}}{\displaystyle \frac{\mathit{\delta }{F}_T}{\mathit{\delta }\mathit{\theta }}},\end{array}$$

where C_a, H_a and ρ_a are the specific heat, height and density of the atmospheric layer, respectively; i is the index representing either land or ocean grid cells; S_a is the absorbed shortwave radiation; I_a is the absorbed longwave radiation; K is the exchange of sensible heat with the surface; L_vP_r is the latent heat released during precipitation, where L_v is the latent heat of vaporization and P_r is precipitation; and L_fP_s is the latent heat released during snowfall, where L_f is the latent heat of fusion and P_s is snowfall. Finally, the last term on the right side of the equation represents the horizontal diffusion of temperature, where R_e is the Earth's radius, θ is the latitude and F_T is the northward flux of dry static energy in the atmospheric layer. All terms in Eq. (1) are in W m⁻².

2.1.1 Radiative fluxes

The model is forced by diurnally averaged solar insolation at the top of the atmosphere ($S_{\mathrm{TOA}}^{\downarrow }$) using the orbital parameter solution from Laskar et al. (2004). The seasonal cycle is driven exclusively by changes in insolation. The absence of a seasonal insolation cycle results in a markedly colder climate. As noted by Bintanja (1997), employing an annual-mean version of their EBM results in insolation no longer being concentrated in the summer months, when lower zenith angles and reduced snow cover enhance the absorption of shortwave radiation. The amount of solar radiation that is reflected at the top of the atmosphere ($S_{\mathrm{TOA}}^{\uparrow }$), transmitted to the surface ($S_{\mathrm{BOA}}^{\downarrow }$) and reflected at the surface ($S_{\mathrm{BOA}}^{\uparrow }$) is calculated using the parameterization from Bintanja (1996). The shortwave parameterization takes into account several atmospheric properties including surface air temperature, surface albedo, solar zenith angles, cloud optical depth and surface height. The shortwave fluxes are computed for both clear-sky and overcast conditions, with the total radiative flux for a given grid cell determined as the weighted average using prescribed cloud cover fractions. In this study, cloud cover fractions over land and ocean are taken from pre-industrial simulations of the Norwegian Earth System Model version 2 (NorESM2; Seland et al., 2020) (Fig. 1a). Daytime-mean solar zenith angles are calculated with equations provided by Balmes and Fu (2020). Of particular importance to the shortwave radiative fluxes is the cloud optical depth parameter (τ). Following Bintanja (1996, 1997) and Stap et al. (2014), τ is kept fixed to a globally and seasonally invariant value. The amount of shortwave radiation that is absorbed by the atmosphere (S_a) and the surface (S_(i)) – either land or ocean – is described below:

$$\begin{array}{}\text{(2)}& {\displaystyle }{S}_{\left(i\right)}=S_{\mathrm{BOA}\left(i\right)}^{\downarrow }-S_{\mathrm{BOA}\left(i\right)}^{\uparrow },\text{(3)}& {\displaystyle }{S}_{\mathrm{a}\left(i\right)}=(S_{\mathrm{TOA}\left(i\right)}^{\downarrow }-S_{\mathrm{TOA}\left(i\right)}^{\uparrow })-S_{\left(i\right)}.\end{array}$$

The outgoing longwave radiative fluxes at the TOA ($I_{\mathrm{TOA}}^{\uparrow }$) and the surface ($I_{\mathrm{BOA}}^{\uparrow }$), together with the incoming longwave flux at the surface ($I_{\mathrm{BOA}}^{\downarrow }$), are also calculated using a radiation parameterization from Bintanja (1996). The longwave parameterization is made a function of surface temperature, surface air temperature, cloud emissivity, surface elevation and atmospheric CO₂ concentrations. As with the shortwave parameterization, the longwave radiative fluxes are calculated separately for clear-sky and overcast conditions. Alterations to the longwave parameterization from Stap et al. (2014), to both increase the climate sensitivity per CO₂ doubling and account for the effects of non-CO₂ greenhouse gases, are maintained. The absorbed longwave radiation fluxes at the surface (I_(i)) and by the atmosphere (I_a) are shown below:

$$\begin{array}{}\text{(4)}& {\displaystyle }{I}_{\left(i\right)}=I_{\mathrm{BOA}\left(i\right)}^{\downarrow }-I_{\mathrm{BOA}\left(i\right)}^{\uparrow },\text{(5)}& {\displaystyle }{I}_{\mathrm{a}\left(i\right)}=I_{\mathrm{BOA}\left(i\right)}^{\uparrow }-(I_{\mathrm{BOA}\left(i\right)}^{\downarrow }+I_{\mathrm{TOA}\left(i\right)}^{\uparrow }).\end{array}$$

Figure 1 (a) Cloud cover fractions over land and ocean taken from pre-industrial simulations of NorESM2 (Seland et al., 2020) for calculations of the shortwave and longwave radiative fluxes. (b) The weighting function from Siler et al. (2018) which determines the fraction of atmospheric heat transport carried out by the Hadley cell. (c) The prescribed vertical ocean velocities for driving ocean circulation, with upwelling from 50° S to 60° N (with an average upwelling rate of 4 m yr⁻¹) and downwelling from 70 to 50° S and from 60 to 80° N. (d) The resulting northward horizontal velocities in the uppermost ocean layer (solid black line), together with the mass transport (dotted grey line), in sverdrups (1×10⁶ m³ s⁻¹).

2.1.2 Turbulent heat fluxes

The aerodynamic bulk relationships are employed to compute the fluxes of sensible heat (K) and evaporation (E) across the atmosphere–ocean and atmosphere–land interfaces:

$$\begin{array}{}\text{(6)}& {\displaystyle }{K}_{\left(i\right)}={\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{a}}{\mathit{\kappa }}_{\left(i\right)}(T_{\left(i\right)}-T_{\mathrm{a}\left(i\right)}),\text{(7)}& {\displaystyle }{E}_{\left(i\right)}={\mathit{\rho }}_{\mathrm{a}}{W}_{\left(i\right)}{\mathit{\kappa }}_{\left(i\right)}\left(Q_{\mathrm{sat}}\right(T_{\left(i\right)})-Q_{\mathrm{a}\left(i\right)}),\end{array}$$

where κ is the turbulent exchange coefficient; T_(i) is the temperature of the surface – either land or ocean; W denotes the surface water availability; Q_sat(T_(i)) is the saturation specific humidity of the surface (as determined by the Clausius–Clapeyron relation); and Q_a is the specific humidity of the overlying atmospheric layer. The influence of wind speed and surface roughness on the turbulent heat exchange is not incorporated into κ, which remains constant across latitudes and seasons. Following Bintanja (1997), W is set to 0.7 and 1.0 over land and ocean, respectively, to reflect reduced water availability over land. The latent heat flux associated with evaporation is simply L_vE.

2.1.3 Hydrological cycle

An atmospheric moisture budget is introduced to parameterize the hydrological cycle and simulate precipitation and snowfall (Fanning and Weaver, 1996; Robinson et al., 2010; Ritz et al., 2011):

$$\begin{array}{}\text{(8)}& {\mathit{\rho }}_{\mathrm{a}}{H}_{\mathrm{a}}{\displaystyle \frac{\mathit{\delta }{Q}_{\mathrm{a}\left(i\right)}}{\mathit{\delta }t}}=(E_{\left(i\right)}-P_{\mathrm{r}\left(i\right)})+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}^{\mathrm{2}}\cos \left(\mathit{\theta }\right)}}{\displaystyle \frac{\mathit{\delta }{F}_Q}{L_{\mathrm{v}}\mathit{\delta }\mathit{\theta }}},\end{array}$$

where the first term on the right side of the equation represents the sources (evaporation) and sinks (precipitation) and the second term represents the horizontal transport of water vapour, in which F_Q is the latent heat associated with the northward moisture flux. Precipitation occurs once the relative humidity, r, exceeds a maximum threshold, r_max:

$$\begin{array}{}\text{(9)}& P_{\mathrm{r}\left(i\right)}=\left\{\begin{array}{ll}\frac{{\mathit{\rho }}_{\mathrm{a}}{H}_{\mathrm{a}}}{\mathrm{3}\cdot \mathrm{\Delta }t}\left(Q_{\mathrm{a}\left(i\right)}-r_{\max }\cdot Q_{\mathrm{sat}}\left(T_{\mathrm{a}\left(i\right)}\right)\right),& r>r_{\max },\\ \mathrm{0},& \text{otherwise},\end{array}\right.\end{array}$$

where Q_sat(T_a) is the saturation specific humidity of the atmospheric layer and Δt is the seconds in 1 d. When the relative humidity exceeds r_max, it leads to precipitation. The turnover time for excess humidity in the atmosphere is set to 3 d. Shorter turnover times produce very large and sporadic contributions of latent heat to the atmospheric column. The amount of snowfall is determined using

$$\begin{array}{}\text{(10)}& P_{\mathrm{s}\left(i\right)}=P_{\mathrm{r}\left(i\right)}\cdot f_{\mathrm{sf}},\end{array}$$

where f_sf is the fraction of precipitation that falls as snow, which is parameterized as a function of surface air temperature (Harvey, 1988):

$$\begin{array}{}\text{(11)}& f_{\mathrm{sf}\left(i\right)}=\left\{\begin{array}{ll}\mathrm{1},& T_{\mathrm{a}\left(i\right)}^{\prime }<\mathrm{260}\mathrm{K},\\ \mathrm{0.05}(\mathrm{280}-T_{\mathrm{a}\left(i\right)}^{\prime }),& \mathrm{260}\mathrm{K}\le T_{\mathrm{a}\left(i\right)}^{\prime }\le \mathrm{280}\mathrm{K},\\ \mathrm{0},& T_{\mathrm{a}\left(i\right)}^{\prime }>\mathrm{280}\mathrm{K},\end{array}\right.\end{array}$$

where $T_{\mathrm{a}}^{\prime }$ is the surface air temperature corrected for the zonal-mean elevation with a global mean lapse rate of −6.5 K km⁻¹. In the current version of the model, this zonal-mean elevation is prescribed (see Sect. 2.4 and Table 2). In the future, we intend to make the zonal-mean elevation dependent on a coupled ice sheet model. The expression for f_sf is taken from Harvey (1988) as the fractional area of a grid box over which precipitation falls as snow, based on meteorological station data. Therefore, rather than assuming a uniform distribution of snowfall across each grid box, this parameterization allows for only a portion of the land or ocean surface to be snow-covered. As precipitation is assumed to fall uniformly over each grid box, however, this geographic fraction also represents the overall proportion of precipitation that is converted into snow.

2.1.4 Atmospheric heat transport

The division of the model into land and ocean raises the question of how to parameterize the “zonal mixing” of air belonging to the same latitudes but overlying each surface type. We adopt an “infinite wind” mixing scenario (Peng et al., 1987; Bintanja, 1997), in which atmospheric temperatures and humidities over land and ocean are both set equal to the zonal mean ($\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ and $\stackrel{\mathrm{‾}}{Q_{\mathrm{a}}}$, respectively) at the end of each model time step, according to the land fraction at each grid cell. The assumption of infinite wind mixing between atmosphere over land and ocean is based on previous studies that find that a more realistic representation of zonal mixing produces small differences in model sensitivity (Harvey and Schneider, 1985; Thompson and Schneider, 1979). In other words, one atmospheric layer effectively overlies the entire surface. Following the zonal mixing of temperature and humidity, meridional atmospheric heat transport is parameterized as a diffusive process along horizontal gradients in moist static energy ($\stackrel{\mathrm{‾}}{m}$), expressed as $\stackrel{\mathrm{‾}}{m}=c_{\mathrm{a}}\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}+L_{\mathrm{v}}\stackrel{\mathrm{‾}}{Q_{\mathrm{a}}}$, where $c_{\mathrm{a}}\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ and $L_{\mathrm{v}}\stackrel{\mathrm{‾}}{Q_{\mathrm{a}}}$ represent the dry static and moist components of atmospheric heat content, respectively. Moist static EBMs behave differently to classic (dry) EBMs (Hwang and Frierson, 2010; Feldl and Merlis, 2021) and have shown success at emulating the response of more comprehensive GCMs to climate forcings (Hwang and Frierson, 2010; Roe et al., 2015). Consequently, the total northward transport of energy (F_total) within the atmospheric layer is described as follows:

$$\begin{array}{}\text{(12)}& F_{\mathrm{total}}=\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{a}}{H}_{\mathrm{a}}{D}_{\mathrm{a}}{\displaystyle \frac{\mathit{\delta }\stackrel{\mathrm{‾}}{m}}{R_{\mathrm{e}}\mathit{\delta }\mathit{\theta }}},\end{array}$$

where D_a is the atmospheric diffusion coefficient. The total atmospheric heat transport is divided into the dry static (F_T) and latent component (F_Q). F_T contributes directly to the heating of the atmospheric layer, whereas F_Q transports water vapour within the simplified hydrological cycle (Eq. 8), which is translated into heating of the atmospheric column once precipitation occurs. To capture the equatorward transport of latent heat in the tropics, against meridional gradients in $\stackrel{\mathrm{‾}}{m}$, the Hadley cell parameterization introduced by Siler et al. (2018) is included to obtain a more realistic representation of the hydrological cycle. In this parameterization, F_total is partitioned into a Hadley cell (F_HC) and an eddy component (F_eddy) (Siler et al., 2018):

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{F}_{\mathrm{total}}& =F_{\mathrm{HC}}+F_{\mathrm{eddy}}\\ & =F_{\mathrm{total}}\mathrm{\Omega }+F_{\mathrm{total}}[\mathrm{1}-\mathrm{\Omega }],\end{array}\end{array}$$

where Ω is a weighting function (shown in Fig. 1b) that ensures the Hadley cell dominates heat transport in the tropics, whereas eddies control the poleward flux of $\stackrel{\mathrm{‾}}{m}$ in the mid-latitudes and polar regions. Where Ω=0, eddies account for all heat transport via the down-gradient diffusion of both dry static (F_{T_eddy}) and latent (F_{Q_eddy}) heat:

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{F}_{T\text{\_eddy}}=\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{a}}{H}_{\mathrm{a}}{c}_{\mathrm{a}}{D}_{\mathrm{a}}{\displaystyle \frac{\mathit{\delta }\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}}{R_{\mathrm{e}}\mathit{\delta }\mathit{\theta }}}[\mathrm{1}-\mathrm{\Omega }],\text{(15)}& {\displaystyle }& {\displaystyle }{F}_{Q\text{\_eddy}}=\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{a}}{H}_{\mathrm{a}}{L}_{\mathrm{v}}{D}_{\mathrm{a}}{\displaystyle \frac{\mathit{\delta }\stackrel{\mathrm{‾}}{Q_{\mathrm{a}}}}{R_{\mathrm{e}}\mathit{\delta }\mathit{\theta }}}[\mathrm{1}-\mathrm{\Omega }].\end{array}$$

In the Hadley cell, the poleward flux of $\stackrel{\mathrm{‾}}{m}$ in its upper branch slightly exceeds the equatorward flux in its lower branch (Hartmann, 2015). Noting that meridional gradients in $\stackrel{\mathrm{‾}}{m}$ are relatively flat within this upper branch, Siler et al. (2018) approximate this difference in $\stackrel{\mathrm{‾}}{m}$ between the upper and lower branch as $g=\mathit{\lambda }\cdot {\stackrel{\mathrm{‾}}{m}}_{\mathrm{eq}}-\stackrel{\mathrm{‾}}{m}$, where ${\stackrel{\mathrm{‾}}{m}}_{\mathrm{eq}}$ is the near-surface moist static energy at the Equator and λ is the fractional increase in moist static energy in the upper branch of the Hadley cell relative to ${\stackrel{\mathrm{‾}}{m}}_{\mathrm{eq}}$. In Siler et al. (2018) λ is set to 1.06, but for this study λ = 1.03 to improve our simulation of precipitation for the pre-industrial period. Consequently, the net poleward transport of heat within the Hadley cell is given as F_HC=g ψ, where ψ is the mass transport within either the upper or the lower branch of the Hadley cell, which can be solved using Eq. (13) as $\mathit{\psi }=\left(F_{\mathrm{total}}\mathrm{\Omega }\right)/g$. Assuming latent heat transport is confined to the lower branch, the dry static (F_{T_HC}) and latent (F_{Q_HC}) contributions to heat transport within the Hadley cell are parameterized using (Siler et al., 2018)

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{F}_{T\text{\_HC}}=\mathit{\psi }[g+L_{\mathrm{v}}{Q}_{\mathrm{a}}],\text{(17)}& {\displaystyle }& {\displaystyle }{F}_{Q\text{\_HC}}=-\mathit{\psi }{L}_{\mathrm{v}}{Q}_{\mathrm{a}}.\end{array}$$

Finally, the Hadley cell and eddy contributions to the total northward transport of dry static (F_T) and latent (F_Q) heat are simply given by

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}{F}_T& =F_{T\text{\_HC}}+F_{T\text{\_eddy}}\\ & =\mathit{\psi }[g+L_{\mathrm{v}}{Q}_{\mathrm{a}}]+\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{a}}{H}_{\mathrm{a}}{c}_{{}_a}{D}_{\mathrm{a}}{\displaystyle \frac{\mathit{\delta }\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}}{R_{\mathrm{e}}\mathit{\delta }\mathit{\theta }}}\\ & \times [\mathrm{1}-\mathrm{\Omega }],\end{array}\text{(19)}& {\displaystyle }\begin{array}{rl}{F}_Q& =F_{Q\text{\_HC}}+F_{Q\text{\_eddy}}\\ & =-\mathit{\psi }{L}_{\mathrm{v}}{Q}_{\mathrm{a}}+\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{a}}{H}_{\mathrm{a}}{L}_{\mathrm{v}}{D}_{\mathrm{a}}{\displaystyle \frac{\mathit{\delta }\stackrel{\mathrm{‾}}{Q_{\mathrm{a}}}}{R_{\mathrm{e}}\mathit{\delta }\mathit{\theta }}}\\ & \times [\mathrm{1}-\mathrm{\Omega }],\end{array}\end{array}$$

where the divergence of F_T and F_Q (converted into a moisture flux by dividing by L_v) is shown on the right side of Eqs. (1) and (8), respectively.

2.2 Land

The energy balance of the land surface is expressed by

$$\begin{array}{}\text{(20)}& C_{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l}}{H}_{\mathrm{l}}{\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{l}}}{\mathit{\delta }t}}=S_{\mathrm{l}}+I_{\mathrm{l}}-K_{\mathrm{l}}-L_{\mathrm{v}}{E}_{\mathrm{l}},\end{array}$$

where C_l, H_l and ρ_l are the specific heat, depth and density of the ground layer, respectively; T_l is the land surface temperature, S_l is the absorbed shortwave radiation; I_l is the absorbed longwave radiation; K_l is the upward sensible heat flux; and E_l is the upward evaporation flux over land.

2.2.1 Snow mass budget

The model evaluates both the proportion of the land surface covered by the snow (f_sc) and the average thickness of the snowpack (d_sc). A fractional area, f_sf, of the precipitation accumulates as snow at the surface according to the hydrological cycle (Sect. 2.1.3). Similarly to the EBM from Harvey (1988), when f_sf is less than the existing area of the snowpack (i.e. f_sf<f_sc), this snowfall gain is redistributed over the larger snow-covered area. Conversely, if f_sf>f_sc, the snow-covered area is expanded to match the area of snowfall (f_sf=f_sc) with the average thickness of the snowpack, d_sc, adjusted to conserve the total mass of snow in each grid cell.

If land surface temperatures (T_l) exceed the melting point of snow (T_K), T_l is reset to T_K and the excess energy is used to melt the snowpack. This rate of melting ($\dot{M}$) is described by

$$\begin{array}{}\text{(21)}& \dot{M}=\left\{\begin{array}{ll}\mathrm{0},& T_{\mathrm{l}}\le T_K,\\ (T_{\mathrm{l}}-T_K)\cdot \frac{C_{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l}}{H}_{\mathrm{l}}}{L_{\mathrm{f}}{\mathit{\rho }}_{\mathrm{ice}}},& T_{\mathrm{l}}>T_K,\end{array}\right.\end{array}$$

where ρ_ice represents the density of ice. Surface melting of the snowpack is evenly distributed between reducing the average thickness, d_sc, and the fractional area of the snowpack, f_sc. Should $\dot{M}$ exceed d_sc, excess melt is reconverted into heating of the land surface.

2.2.2 Land surface albedo

The albedo of the land is determined by the fraction of the surface that is covered in snow, together with the albedo of the snow-covered and snow-free surfaces. The fractional snow cover, f_sc, is assumed to be evenly distributed between two surface types: bare ground and land ice. The snow albedo (α_s) is expressed as a linear function of surface temperature following Bintanja (1997).

$$\begin{array}{}\text{(22)}& {\mathit{\alpha }}_{\mathrm{s}}=\left\{\begin{array}{ll}{\mathit{\alpha }}_{\mathrm{cs}},& T_{\mathrm{l}}^{\prime }<\mathrm{263}\mathrm{K},\\ {\mathit{\alpha }}_{\mathrm{cs}}+[{\mathit{\alpha }}_{\mathrm{ws}}-{\mathit{\alpha }}_{\mathrm{cs}}]\frac{(T_{\mathrm{l}}^{\prime }-\mathrm{263})}{\mathrm{10}},& \mathrm{263}\mathrm{K}\le T_{\mathrm{l}}^{\prime }\le \mathrm{273}\mathrm{K},\\ {\mathit{\alpha }}_{\mathrm{ws}},& T_{\mathrm{l}}^{\prime }>\mathrm{273}\mathrm{K},\end{array}\right.\end{array}$$

where α_cs is the maximum (or “cold”) snow albedo, α_ws is the minimum (or “warm”) snow albedo and $T_{\mathrm{l}}^{\prime }$ is the land surface temperature corrected for zonal-mean elevation. The average albedo over bare ground (α_g) and ice (α_i) is then calculated as the weighted average of the snow-covered and snow-free region.

$$\begin{array}{}\text{(23)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{g}}=f_{\mathrm{sc}}{\mathit{\alpha }}_{\mathrm{s}}+(\mathrm{1}-f_{\mathrm{sc}}){\mathit{\alpha }}_{\mathrm{bg}},\text{(24)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{i}}=f_{\mathrm{sc}}{\mathit{\alpha }}_{\mathrm{s}}+(\mathrm{1}-f_{\mathrm{sc}}){\mathit{\alpha }}_{\mathrm{bi}},\end{array}$$

where α_bg and α_bi are the albedo of bare ground and the albedo of ice without snow cover, respectively. The assumption of a uniform α_bg albedo overlooks the important influence that different vegetation types have on land albedo. In contrast, Bintanja (1997) divides “ice-free” land into present-day distributions of grass and forest cover, though these proportions are held constant over time. Thus, in both approaches, these potentially significant vegetation feedbacks are excluded from Quaternary climate simulations. While including present-day vegetation distribution could improve pre-industrial simulations of ZEMBA, we see limited added value in doing so for studies of orbitally driven climate change. Nonetheless, we recognize that these simplifications in land albedo may affect the strength of albedo feedbacks over land, which could be explored in future applications of the model. The average albedo over land (α_l) is the weighted average of α_g and α_i, depending on the fractional area over which land is covered by ice (f_i):

$$\begin{array}{}\text{(25)}& {\mathit{\alpha }}_{\mathrm{l}}=f_{\mathrm{i}}{\mathit{\alpha }}_{\mathrm{i}}+(\mathrm{1}-f_{\mathrm{i}}){\mathit{\alpha }}_{\mathrm{g}}.\end{array}$$

Additionally, the impact of the solar zenith angle, θ_z, on albedo is accounted for by increasing the land albedo for zenith angles greater than 60° using the following parameterization (Lefebre et al., 2003):

$$\begin{array}{}\text{(26)}& {\mathit{\alpha }}_{\mathrm{l}}=\left\{\begin{array}{ll}{\mathit{\alpha }}_{\mathrm{l}}+max\left\{\mathrm{0};\mathrm{0.32}\cdot \frac{\mathrm{1}}{\mathrm{2}}\left[\frac{\mathrm{3}}{\mathrm{1}+\mathrm{4}\cos {\mathit{\theta }}_{\mathrm{z}}}-\mathrm{1}\right]\right\},& {\mathit{\theta }}_{\mathrm{z}}\ge \mathrm{60}\mathit{°},\\ {\mathit{\alpha }}_{\mathrm{l}},& {\mathit{\theta }}_{\mathrm{z}}<\mathrm{60}\mathit{°},\end{array}\right.\end{array}$$

with the constraint that ${\mathit{\theta }}_{\mathrm{z}}>\mathrm{80}\mathit{°}=\mathrm{80}\mathit{°}$. This adjustment increases land albedo for very high zenith angles in the polar regions.

2.3 Ocean

Meridional heat transport by the oceans is represented by the zonally averaged ocean circulation model from Bintanja (1997). Extending from 70° S to 80° N, the ocean model comprises six layers of increasing thickness with depth (with a total depth of 4000 m). To induce thermohaline circulation, prescribed vertical velocities (shown in Fig. 1c) are used to produce regions of upwelling from 50° S to 60° N and downwelling in the polar regions. These vertical velocities remain constant with depth and drive the poleward flow of water in the uppermost model layer (Fig. 1c) and equatorward flow in the bottom layer, thereby generating a conveyor-belt system of ocean heat transport in each hemisphere. The central point of ocean circulation is placed slightly south of the Equator at 5° S to improve temperature estimates in the north polar regions. The model effectively consists of two overturning cells, with an average upwelling rate of 4 m yr⁻¹ in the upwelling regions (from 50 to 5° S in the southern cell and from 5° S to 60° N in the northern cell). The transport of ocean heat via eddies and gyres is represented as a diffusive process, along horizontal gradients in surface ocean temperature. Outside of the ocean circulation basin, where ocean fractions are greater than zero (e.g. from 80° S to 90° N), a passive mixed layer of 100 m depth exchanges radiative and turbulent heat with the atmosphere. The temporal evolution of ocean temperature, T_o, is described by (Bintanja, 1997)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{C}_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}{H}_{\mathrm{o}}{\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{o}}}{\mathit{\delta }t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}^{\mathrm{2}}f\cos \mathit{\theta }}}{\displaystyle \frac{\mathit{\delta }{F}_{\mathrm{ov}}}{\mathit{\delta }\mathit{\theta }}}+C_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}{H}_{\mathrm{o}}{\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }z}}\left(w{T}_{\mathrm{o}}\right)\\ {\displaystyle }& {\displaystyle }=[S_{\mathrm{o}}+I_{\mathrm{o}}-K_{\mathrm{o}}-L_{\mathrm{v}}{E}_{\mathrm{o}}]+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}^{\mathrm{2}}f\cos \mathit{\theta }}}{\displaystyle \frac{F_{\mathrm{eg}}}{\mathit{\delta }\mathit{\theta }}}+\\ \text{(27)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}^{\mathrm{2}}f\cos \mathit{\theta }}}{\displaystyle \frac{\mathit{\delta }{F}_{\mathrm{i}}}{\mathit{\delta }\mathit{\theta }}}+C_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}{H}_{\mathrm{o}}{\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }z}}\left(D_z{\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{o}}}{\mathit{\delta }z}}\right),\end{array}$$

where C_o, ρ_o and H_o are the density, specific heat and depth of each ocean layer; z is the vertical coordinate; f is the fractional width of the ocean basin; F_ov is the northward advective heat flux in the top and bottom ocean layer due to overturning; w is the prescribed vertical velocity; S_o is absorbed shortwave radiation; I_o is absorbed longwave radiative flux; K_o is the exchange of sensible heat with the atmosphere; E_o is evaporation at the ocean surface; F_eg is the northward diffusive heat flux at the ocean surface due to eddies and gyres; F_i is the northward diffusive heat flux in the ocean interior; and D_z is the coefficient for vertical heat diffusion. The second term on the left side of the equation represents the divergence of horizontal heat advection which pertains to the top and bottom layer of the ocean model, and the third term represents the divergence of vertical heat advection. On the right side, the first to fourth terms represent the radiative and turbulent heat fluxes and the fifth term is the divergence of horizontal heat diffusion (via eddies and gyres), all of which apply to the uppermost ocean layer. The sixth term on the right side – the divergence of horizontal heat diffusion in the ocean interior – applies to every layer except the surface. Finally, the last term on the right side of the equation represents the divergence of vertical heat diffusion. The northward heat flux associated with F_ov, F_eg and F_i is described by

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{ov}}=\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}f\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{o}}{C}_{\mathrm{o}}{H}_{\mathrm{o}}u{T}_{\mathrm{o}},\text{(29)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{eg}}=-\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}f\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{o}}{C}_{\mathrm{o}}{H}_{\mathrm{o}}{D}_{\mathrm{o}}{\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{o}}}{R_{\mathrm{e}}\mathit{\delta }\mathit{\theta }}},\text{(30)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{i}}=-\mathrm{2}\mathit{\pi }{R}_{\mathrm{e}}f\cos \left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{o}}{C}_{\mathrm{o}}{H}_{\mathrm{o}}{D}_{\mathrm{i}}{\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{o}}}{R_{\mathrm{e}}\mathit{\delta }\mathit{\theta }}},\end{array}$$

where u is the horizontal ocean velocity, D_o is the diffusion coefficient related to eddy and gyre transport at the surface and D_i is the coefficient related to horizontal diffusion in the ocean interior. The distribution of u is calculated using the continuity equation from Bintanja (1997, p. 9) and is shown in Fig. 1d.

2.3.1 Sea ice

The inclusion of sea ice is important for capturing the seasonal range of surface temperatures at the higher latitudes. By reflecting the majority of incoming solar radiation and reducing heat exchange between the atmosphere and ocean surface, the presence of sea ice reduces the effective thermal inertia of the ocean–atmosphere system. Therefore, ZEMBA includes a simple sea ice model to simulate the latitudinal distribution of sea ice and its modification of ocean albedo. The sea ice model does not account for variations in sea ice thickness or sea ice drifting. Sea ice of a prescribed thickness forms or melts when surface ocean temperatures drop below or exceed a critical threshold (T_fo) – the freezing temperature of seawater. The heat flux available for the formation or melting of sea ice, Q_si, is governed by (Gildor and Tziperman, 2001)

$$\begin{array}{}\text{(31)}& Q_{\mathrm{si}}={\displaystyle \frac{C_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}{H}_{\mathrm{o}}{A}_{\mathrm{o}}}{\mathrm{\Delta }t}}(T_{\mathrm{fo}}-T_{\mathrm{o}}),\end{array}$$

where A_o is the surface area of the ocean. The heat flux is converted into changes in sea ice volume, V_si, as follows:

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\mathit{\delta }{V}_{\mathrm{si}}}{\mathit{\delta }t}}={\displaystyle \frac{Q_{\mathrm{si}}}{{\mathit{\rho }}_{\mathrm{ice}}{L}_{\mathrm{f}}}}+{\displaystyle \frac{P_{\mathrm{s}}\cdot A_{\mathrm{si}}}{{\mathit{\rho }}_{\mathrm{ice}}}},\end{array}$$

where the last term on the right side of the equation represents the snowfall contribution to sea ice volume, with A_si being the surface area of the sea ice cover. Sea ice volume is then converted into sea ice areal extent by assuming a constant sea ice thickness (d_si), which is set to 2 m. When sea ice forms, the temperature of the underlying ocean is reset to T_fo, and when sea ice covers the entire ocean surface (i.e. A_si=A_o), the temperature of the sea ice can drop below the freezing point, while the upper ocean layer remains at T_fo. Variations in the temperature of the surface ocean layer (and thereby sea ice) are determined by the surface energy balance, in addition to advective and diffusive ocean heat fluxes within the ocean circulation model.

2.3.2 Ocean albedo

The albedo of the open ocean, α_op, is parameterized as a function of the solar zenith angle, as derived from aircraft observations (Taylor et al., 1996):

$$\begin{array}{}\text{(33)}& {\mathit{\alpha }}_{\mathrm{op}}={\displaystyle \frac{\mathrm{0.037}}{\mathrm{1.1}{\mathit{\theta }}_{\mathrm{z}}^{\mathrm{1.4}}+\mathrm{0.15}}},\end{array}$$

with the average albedo for the ocean surface, α_o, described by

$$\begin{array}{}\text{(34)}& {\mathit{\alpha }}_{\mathrm{o}}={\mathit{\alpha }}_{\mathrm{si}}{f}_{\mathrm{si}}+{\mathit{\alpha }}_{\mathrm{op}}(\mathrm{1}-f_{\mathrm{si}}),\end{array}$$

where α_si is the albedo of the sea ice and f_si is the fraction of the ocean surface covered by sea ice.

2.4 Numerical details

The model has the option to be simulated at a 1, 2.5 and 5° resolution. All equations are solved using the forward Euler method, with the majority of calculations performed using a 1 d time step. The exception is the atmospheric and ocean heat transport processes, which may be solved using a shorter time step for numerical stability reasons depending on the choice of model resolution. For all experiments described in the subsequent sections, a model resolution of 5° is chosen. Cloud emissivity is set to 1. The percentage of land and ocean cover for each zonal band is taken from the ICE-6G_C dataset (Argus et al., 2014; Peltier et al., 2014), although a mixed ocean layer is assumed to cover the entire surface from 80–90°) and land occupies the whole area poleward of 75°. Values for key model parameters are summarized in Table 1. These are based on values used in previous studies using an EBM which formed the basis of ZEMBA (Bintanja, 1997; Stap et al., 2014) but with small adjustments to improve the simulated pre-industrial zonal-mean temperature. The coefficient for atmospheric heat transport (D_a) has been modified in both hemispheres to improve the simulated polar temperature. Additional sensitivity experiments are presented in Appendix A. The choice of pre-industrial cloud cover can have a significant impact on the simulated climate (Appendix A1). Furthermore, sensitivity experiments conducted for all key model parameters reveal that ZEMBA is particularly sensitive to cloud cover parameters (Appendix A2). Notably, ZEMBA shows a strong sensitivity to the globally averaged cloud optical depth (τ), which has been used as a tuning parameter to adjust the radiation budget to match that of the present day (Bintanja, 1997; Stap et al., 2014). A comprehensive list of all model variables, parameters and constants is provided in Appendix B.

Table 1 Selection of important parameters used in the atmospheric, land and ocean components of the model.

3 Results

We present equilibrium simulations of the model for the pre-industrial period (ZEMBA-PI), a doubling of atmospheric CO₂ conditions relative to the pre-industrial period (ZEMBA-2×CO₂), a 2 % increase in solar insolation (ZEMBA-S_o+2 %) and the Last Glacial Maximum (ZEMBA-LGM). The model takes approximately 3000 model years to reach an equilibrium due to the inclusion of the ocean model with a large heat capacity. Key model parameters remain the same between experiments (Table 1), with the exception of the insolation forcing, atmospheric CO₂ concentration and land ice extent (Table 2). The following sections describe output from ZEMBA prior to coupling with an interactive ice sheet model, with fixed land ice fractions and zonal-mean elevations over land taken from ICE-6G_C (Argus et al., 2014; Peltier et al., 2014).

Table 2 Boundary conditions for ZEMBA experiments. PD denotes present day.

3.1 Pre-industrial simulation

For the ZEMBA-PI experiment, the model is forced with insolation using the present-day orbital parameters (Laskar et al., 2004) and an atmospheric CO₂ concentration of 284 ppm. To evaluate the accuracy of the ZEMBA-PI model output, we compare it against a selection of zonally averaged climate variables from a pre-industrial simulation of the Norwegian Earth System Model version 2 (NorESM2) (Seland et al., 2020). In addition, we compare ZEMBA-PI against the ERA5 atmospheric reanalysis product averaged between 1940 and 1970 (Hersbach et al., 2023), which belongs to a period when CO₂ concentrations were ∼ 26–45 ppm higher than the pre-industrial levels but provides one of the earliest observational constraints on the full-field climate. As can be seen in the subsequent figures (Figs. 2–7), zonally averaged climate variables are very similar between NorESM2 and ERA5 1940–1970.

The zonal-mean surface air temperatures are shown in Fig. 2. For annual-mean temperatures (Fig. 2a), the latitudinal structure corresponds nicely between ZEMBA-PI, NorESM2 and ERA 1940–1970. Moreover, the DJF mean (Fig. 2b) and JJA mean (Fig. 2c) demonstrate that the model successfully captures the seasonal amplitude of temperatures at the higher latitudes in accordance with both NorESM2 and ERA5 1940–1970. The globally averaged surface air temperature from ZEMBA-PI is 13.87 °C, compared with 13.78 °C from NorESM2 and 13.72 °C from ERA5 1940–1970 (Table 3). When the ZEMBA-PI experiment is compared only to NorESM2, the difference in annual-mean temperatures never exceeds ∼ 4 °C (Fig. 2d). The same is true for ERA5 1940–1970, with the exception of 70–90° N, as ERA5 produces higher temperatures in this region due to much higher winter temperatures (Fig. 2b), which is presumably in response to higher CO₂ levels. Between 60 and 30° S, ZEMBA-PI appears to consistently overestimate temperatures, both annually (Fig. 2d) and seasonally (Fig. 2e–f), by as much as ∼ 5 °C when compared to NorESM2 and ERA5 1940–1970. Between 30 and 60° N, on the other hand, the annual-mean temperatures agree well amongst the models and reanalysis, but ZEMBA-PI seems to underestimate the seasonal range in temperatures at these latitudes. Overall, despite some biases, ZEMBA captures both the annual mean and the seasonal range of surface air temperatures with good accuracy for the pre-industrial period.

Figure 2 (a–c) The annual (a), December–January–February (b) and June–July–August (c) average of zonal-mean surface air temperature for the pre-industrial period (PI), as simulated by ZEMBA (black lines) in comparison to NorESM2 (blue lines) and the ERA5 climatology from 1940 to 1970 (red lines). (d–f) The difference between ZEMBA and the other models and observations.

Figure 3 shows climate variables related to the hydrological cycle including precipitation, snowfall and evaporation. Given the simplicity of ZEMBA, precipitation rates from the ZEMBA-PI experiment are well captured when compared to NorESM2 and ERA5 1940–1970 (Fig. 3a). The inclusion of the Hadley cell parameterization from Siler et al. (2018) generates a convergence of moisture and strong precipitation rates near the Equator. The largest difference between ZEMBA-PI and the other models and reanalysis is located around the Equator (Fig. 3c), as the precipitation maximum is located at the Equator for ZEMBA-PI, whereas the maxima for NorESM2 and ERA5 1940–1970 are located north of the Equator in accordance with the mean position of the Intertropical Convergence Zone (ITCZ). In the mid-latitudes and polar regions, precipitation rates from ZEMBA-PI are in close agreement (within 1 mm d⁻¹) with those simulated by the complex models. The snowfall rates (Fig. 3) from ZEMBA-PI over the NH correspond nicely to NorESM2 and ERA5 1940–1970. For the Southern Hemisphere (SH), on the other hand, while ZEMBA-PI captures the location of maximum snowfall at 60° S, it appears to underestimate the snowfall rate by half (Fig. 3c). As for evaporation, the zonally and annually averaged fluxes from ZEMBA-PI are very similar to NorESM2 and ERA5 (Fig. 3b), with the differences between them never exceeding 1 mm d⁻¹ (Fig. 3d).

Figure 3 (a–b) The annual- and zonal-mean rates of precipitation and snowfall (a) and evaporation (b), as simulated by ZEMBA (black lines) in comparison to NorESM2 (blue lines) and the ERA5 1940–1970 climatology (red lines). Precipitation is shown in solid lines, and snowfall is shown in dotted lines (a, c). (c–d) The difference between ZEMBA and the other models and reanalysis.

The surface and planetary albedo compares favourably between ZEMBA-PI, NorESM2 and ERA5 1940–1970 (Fig. 4). Most importantly, ZEMBA-PI reproduces the poleward enhancement of albedo at the higher latitudes due to the presence of snow and sea ice cover. The difference in surface and planetary albedo as simulated by ZEMBA-PI in comparison to the other models and reanalysis never exceeds 0.2 (Fig. 4c–d). One limitation of the shortwave parameterization used in ZEMBA is the overestimation of the planetary albedo by ∼ 0.1 at the polar latitudes (Fig. 4c) as noted by Bintanja (1997). However, as the planetary albedo is slightly underestimated at the mid-latitudes, the globally averaged planetary albedo is very similar between ZEMBA (0.31), NorESM2 (0.32) and ERA5 1940–1970 (0.31) (Table 3). At the surface, the onset of higher surface albedo agrees well with NorESM2 and ERA5 1940–1970 in both hemispheres. Overall, the global mean surface albedo for ZEMBA-PI is 0.15, which is in agreement with both NorESM2 and ERA5 1940–1970 (Table 3).

Figure 4 (a–b) The annual- and zonal-mean planetary (a) and surface (b) albedo for the pre-industrial period (PI), as simulated by ZEMBA (black lines) in comparison to NorESM2 (blue lines) and the ERA5 climatology from 1940 to 1970 (red lines). (c–d) The difference between ZEMBA and the other models and observations.

Figure 5 shows the seasonal cycle in the areal extent of snow cover over land and sea ice coverage. Over land in the NH (Fig. 5a), the ZEMBA-PI experiment underestimates the winter maximum in snow coverage by $>\mathrm{1}\times {\mathrm{10}}^{\mathrm{13}}$ m² and thereby simulates a smaller seasonal amplitude in snow cover when compared to NorESM2 and ERA5 1940–1970 reanalysis. As for the SH (Fig. 5a), all models and reanalysis have similar areal extents in snow cover and produce negligible seasonal variations. The ZEMBA-PI experiment shows most success in simulating sea ice cover in the NH (Fig. 5b), capturing both the amplitude and the phase of sea ice cover changes in reference to NorESM2 and ERA5 1940–1970 (Fig. 5b). For sea ice cover in the SH (Fig. 5b), however, the ZEMBA-PI experiment underestimates the seasonal amplitude.

Figure 5 Monthly variations in the areal extent of snow coverage over land (a) and sea ice coverage over the ocean (b) for the pre-industrial period, as simulated by ZEMBA (black lines) in comparison to NorESM2 (blue lines) and the ERA5 1940–1970 climatology (red lines) for the Northern Hemisphere (solid lines) and the Southern Hemisphere (dotted lines). For snow cover over land, monthly averaged ERA5 reanalysis is taken over a shorter period from 1950 to 1970 (Muñoz Sabater, 2019) due to data availability.

The annual- and zonal-mean radiative fluxes exchanged at the TOA and the surface are shown in Fig. 6. For the TOA, the amount of absorbed shortwave radiation (ASR) estimated using the shortwave parameterization from ZEMBA (see Sect. 2.1.1) compares favourably with NorESM2 and ERA5 1940–1970, albeit with a slight underestimation of ASR in the polar latitudes by ∼ 10–20 W m⁻² (Fig. 6a). The largest differences between ZEMBA-PI and the other models and reanalysis reside in the outgoing longwave flux (OLR) at the TOA (Fig. 6a). In particular, there is a pronounced overestimation of OLR in the tropics (of up to ∼ 30 W m⁻²) and a pronounced underestimation of OLR in the polar latitudes, especially in the Northern Hemisphere (of up to 50 W m⁻²). Consequently, these differences in OLR are reflected in the net radiative flux (NET) received at the TOA, with ZEMBA-PI receiving ∼ 25 W m⁻² less NET energy at the Equator and up to ∼ 35 W m⁻² more NET energy in the high northern latitudes (Fig. 6c). As for the surface radiative fluxes (Fig. 6b), the ASR, OLR and NET radiative fluxes are quite similar but tend to be slightly underestimated in ZEMBA-PI when compared to the other datasets, with the difference in the NET fluxes never exceeding ∼ 25 W m⁻² (Fig. 6d). Overall, despite some clear discrepancies in the ASR and OLR, the ZEMBA-PI experiment produces net radiative fluxes at the TOA and the surface which compare favourably with more complex models, generally falling within ∼ 25 W m⁻² of those simulated by NorESM2 and ERA5 (Fig. 6c–d).

Table 3 Selection of annual-mean variables from a pre-industrial simulation of ZEMBA in comparison to NorESM2 and ERA5 1940–1970 climatology.

Figure 6 (a–b) The annual- and zonal-mean radiative fluxes at the top of the atmosphere (a) and the surface (b), as simulated by ZEMBA (black lines) in comparison to NorESM2 (blue lines) and the ERA5 1940–1970 climatology (red lines). Shown are the absorbed shortwave radiation (ASR: solid lines), outgoing longwave radiation (OLR: dotted lines) and net radiation (NET: dash-dotted lines) radiative fluxes. The outgoing longwave radiation at the surface refers to net upward longwave radiation, i.e. the upward longwave radiative flux at the surface minus the downward longwave radiative flux at the bottom of the atmospheric layer. (c–d) The difference in the net radiative flux between ZEMBA and the other models and reanalysis.

The simulated northward heat transport via the atmosphere and ocean is depicted in Fig. 7 in reference to NorESM2. We note that NorESM2 heat transport values replicate those estimated from 2000 to 2014 using ERA-Interim reanalysis (Trenberth and Fasullo, 2017), including total heat transport exceeding 5.5 PW in each hemisphere and ocean heat transport peaking at around 2 PW at 15° N. The total heat transport in ZEMBA-PI is lower than in NorESM2 in each hemisphere (Fig. 7a) because of reduced heat transport in the atmospheric layer (Fig. 7b). Ocean heat transport corresponds nicely to that inferred from NorESM2 (Fig. 7b), with a maximum value of ∼ 1.5 and ∼ 1 PW in the NH and SH, respectively. However, the location of maximum ocean heat transport is located at ∼ 30° N/S in the ZEMBA-PI experiment compared to ∼ 15° N/S for NorESM2. While the maximum atmospheric heat transport is underestimated in ZEMBA-PI, the location of maximum heat transport corresponds nicely to NorESM2 at ∼ 45° N/S. For the ZEMBA-PI experiment, latent heat transport contributes significantly to the poleward flux of atmospheric heat at the mid-latitudes, whereas dry static heat transport dominates the polar regions, in keeping with NorESM2 (Fig. 7c). Moreover, dry static transport peaks at ∼ 15° N/S and latent heat transport peaks at ∼ 40° N/S for both ZEMBA-PI and NorESM2. The inclusion of the Hadley cell parameterization (Siler et al., 2018) successfully produces the equatorward flux of latent heat seen in the tropics. Overall, while dry static and latent heat transport in the atmosphere is underestimated in the ZEMBA-PI experiment compared to NorESM2, the overall configuration of atmospheric heat transport is very similar between the models.

Figure 7 The northward transport of atmospheric and ocean heat in the pre-industrial period, as simulated by ZEMBA (black lines) in comparison to NorESM2 (blue lines). Shown is the total heat transport (a), the atmospheric and ocean components (b), and the partition of atmospheric heat transport into the dry static and latent components (c). For NorESM2, annual-mean heat transport is inferred by assuming the system is in equilibrium and the heat transport is equal to energy imbalance at any latitude. More specifically, the energy imbalance at the top of the atmosphere (for total heat transport), in the atmosphere (for atmospheric heat transport) or at the surface (for ocean heat transport) is integrated from the South Pole.

3.2 2xCO₂ and +2 % insolation

A common method to test the sensitivity of climate models is to impose changes in atmospheric CO₂ concentrations or solar insolation (S₀). Therefore, the temperature response to a doubling of atmospheric CO₂ concentrations (ZEMBA-2×CO₂) and a 2 % increase in the solar constant (ZEMBA-S_o+2 %) is shown in Fig. 8, keeping all other boundary conditions the same as in the ZEMBA-PI experiment. It should be noted that land ice fractions and elevations are kept fixed for these experiments. For a doubling of atmospheric CO₂ concentrations, the global mean temperature is 3.6 °C higher. The most notable feature is that the NH is significantly more sensitive than the SH, reaching an annual-mean rise of >15 °C in the high northern latitudes, compared to <8 °C over Antarctica (Fig. 8a). In addition, temperature changes are strongest during the winter months of both hemispheres (Fig. 8b). The response to a 2 % increase in the solar constant is very similar to a doubling of atmospheric CO₂ concentrations, with an equivalent rise in global mean surface air temperature of 3.3 °C but with more muted warming over Antarctica. In reference to other works, comparisons are made complicated by the fact that GCM simulations involving a doubling or quadrupling of atmospheric CO₂ are often not run to their equilibrium climate state due to considerable computing times. However, the global mean warming from the ZEMBA-2×CO₂ experiment (3.6 °C), known as the equilibrium climate sensitivity (ECS), fits comfortably within the “likely” range of 1.5–4.5 °C estimated in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Collins et al., 2013). In addition, polar amplification in surface warming that is strongest in the NH and during the winter months is consistent with both historical observations (1979–2014) and GCMs responding to an abrupt quadrupling of the atmospheric CO₂ concentration (Hahn et al., 2021).

Figure 8 Changes in surface air temperature (relative to the PI period) following a doubling of the atmospheric CO₂ concentration from 284 to 568 ppm (a, b; 2xCO₂) and a 2 % increase in the solar constant (c, d; S_o+2 %). (a, c) Changes in the annual-mean temperature, which are also shown normalized by the global mean warming (on the right-hand side). (b, d) Changes in the seasonal cycle in temperature.

3.3 Last Glacial Maximum

To test the ability of ZEMBA to simulate climates other than the pre-industrial period, we perform a simulation of the LGM (hereafter the ZEMBA-LGM experiment). The LGM remains a natural focus for climate models constraining their climate sensitivity because it was the most recent cold extreme when atmospheric CO₂ concentrations were ∼ 100 ppm lower than in the pre-industrial period (Bereiter et al., 2015) and continental ice sheets reached their maximum extent over North America and Fennoscandia. For LGM boundary conditions, the model is forced with insolation using the orbital parameters from 21 ka (Laskar et al., 2004), prescribed changes in land elevation and land ice fractions from ICE-6G_C (Argus et al., 2014; Peltier et al., 2014), and an atmospheric CO₂ concentration of 184 ppm (Bereiter et al., 2015) (see Table 2). The differences between the LGM and PI climates, as simulated by ZEMBA, are compared to an ensemble of state-of-the-art climate models which contributed to PMIP3 and PMIP4 compiled by Kageyama et al. (2021). In addition, we compare our simulation to a recent LGM reconstruction from Annan et al. (2022) which combines an ensemble of climate model simulations with proxy-based estimates of surface temperature using a data assimilation approach. In the ZEMBA-LGM experiment shown in Fig. 9, there are no changes to the strength or the configuration of ocean circulation.

Figure 9 (a, b) The LGM minus PI (LGM − PI) surface air temperatures (a) and precipitation (b) as simulated by ZEMBA (solid black line) in comparison to the multi-model mean from PMIP3 (dash-dotted orange line) and PMIP4 (dash-dotted yellow line) and (in panel a only) the data assimilation products from Annan et al. (2022) (dash-dotted blue line). The yellow- and blue-shaded areas represent the range of LGM − PI surface air temperatures and precipitation as reconstructed by PMIP3–PMIP4 and Annan et al. (2022), respectively. (c, d) The differences in northward heat transport between the LGM and PI experiments for ZEMBA, including total heat transport (solid black line), atmospheric heat transport (solid red line), ocean heat transport (solid blue line), the decomposition of ocean heat transport into ocean overturning (dotted blue line) and eddy and gyre transport (dashed blue line) (c) and the decomposition of atmospheric heat transport (solid red line) into dry static (dotted red line) and latent (dashed green line) heat transport (d).

The changes in annual- and zonal-mean surface air temperatures for the ZEMBA-LGM simulation are shown in Fig. 9a. The temperature decrease around the tropics (30° S–30° N) of −1.95 °C is comparable to the other reconstructions but falls slightly on the lower end of estimated cooling (Table 4). In the extratropics, however, ZEMBA appears to generate cooling which is slightly too strong in the NH and too mild in the SH when compared to the multi-model averages from PMIP3, PMIP4 and the Annan et al. (2022) reconstruction. The latitudinal distribution of cooling always falls within the range simulated across the PMIP3–PMIP4 model ensemble (yellow-shaded area, Fig. 9), but this is not always the case at the higher latitudes for the Annan et al. (2022) reconstruction (blue-shaded area, Fig. 9), which is constrained by proxies for surface temperatures. Overall, despite some discrepancies in the higher latitudes, the global mean cooling from the ZEMBA-LGM experiment is −4.11 °C compared to −4.71 °C from PMIP3, −4.77 °C from PMIP4 and −4.46 °C from Annan et al. (2022) (Table 4). For LGM − PI precipitation rates (Fig. 9b), ZEMBA captures the widespread decrease in precipitation, which again falls within the range estimated by the PMIP3–PMIP4 model ensemble. The reductions are largest in the NH extratropics, similarly to surface air temperatures. However, there are a number of discrepancies between ZEMBA LGM − PI precipitation rates and the PMIP3 and PMIP4 model ensembles, such as at 5–10° N, where the PMIP ensemble generates a much stronger reduction in precipitation, and at 30–10° S, where the PMIP ensemble simulates higher zonal-mean precipitation in the LGM due to increased precipitation in the subtropical Pacific Ocean (Kageyama et al., 2021). Overall, the global mean decrease in precipitation is −0.24 mm d⁻¹ for ZEMBA compared to −0.29 mm d⁻¹ in PMIP3 and −0.33 mm d⁻¹ in PMIP4 (Table 5).

The differences in northward heat transport between the ZEMBA-LGM and ZEMBA-PI experiments are shown in Fig. 9c–d. Starting with total heat transport (Fig. 9c), the most notable feature is a significant peak in northward heat transport (of ∼ 0.5 PW) at 45° N, followed by a significant trough (of ∼ 0.3 PW) at 60° N. The peak is associated with greater atmospheric heat transport in the LGM experiment (Fig. 9c–d), which is almost entirely because of larger fluxes of dry static energy in the atmospheric layer (Fig. 9d). The reduced northward heat transport at 60° N, on the other hand, is due to a significant 0.5 PW decrease in ocean heat transport (Fig. 9c), associated with sea ice expansion in the LGM experiment. As the surface ocean layer underlying sea ice rests at the freezing point of seawater (see Sect. 2.3.1), there is effectively zero meridional ocean heat fluxes at latitudes covered by sea ice. Consequently, sea ice expansion is accompanied by large drops in ocean heat transport, via both overturning and eddy and gyre transport, at latitudes which were previously free of sea ice in the PI experiment.

Figure 10 The LGM − PI surface air temperatures (a) and precipitation rates (b) for a series of experiments, which involve progressively adding LGM land ice fractions and zonal-mean land elevations (dash-dotted blue line), LGM CO₂ concentrations (dash-dotted red line), LGM insolation (black line), a shift in the mid-point of ocean circulation from 5 to 15° S (dash-dotted black line), a 75 % decrease in ocean overturning (dash-dotted green line) strength, and a 50 % decrease in ocean overturning strength (dash-dotted purple line).

In addition to the standard ZEMBA-LGM simulation, we recreate the sensitivity experiment from Bintanja and Oerlemans (1996) by progressively adding each of the LGM boundary conditions (shown in Table 2) followed by changes to the strength of ocean circulation (Fig. 10). Using the EBM from Bintanja (1997), Stap et al. (2014) also found that introducing a mechanism to shift the mid-point of ocean circulation further south during glacial conditions was necessary to produce surface temperatures in closer agreement with observation. Therefore, we additionally investigate the impact of changing the mid-point of ocean circulation from 5 to 15° S. To summarize, we progressively add each of the following:

• ice – LGM land ice fractions and zonal-mean land elevations;

• CO₂ – LGM CO₂ concentrations;

• Inso – LGM insolation forcing, which is the same as in the standard ZEMBA-LGM experiment;

• Oc: 15° S – a shift in the mid-point of ocean circulation from 5 to 15° S;

• 75 % Ov – a 25 % reduction in the strength of ocean overturning;

• 50 % Ov – another 25 % (50 % total) reduction in the strength of ocean overturning.

The changes in temperature and precipitation – averaged over the tropics, extratropics and globally – following the addition of each boundary condition are shown in Tables 4 and 5, respectively. Firstly, the addition of LGM land ice fractions and elevation causes a large decrease in temperature (−5.06 °C) and precipitation (−0.22 mm d⁻¹) in the NH extratropics, whereas the effects in the tropics and the SH extratropics are minimal. A subsequent reduction in atmospheric CO₂ concentrations generates widespread cooling, drying and a further reduction in global mean temperature of −2.08 °C, the latter of which is stronger than the initial cooling caused by the LGM ice sheets (−1.69 °C). The subsequent addition of the LGM insolation forcing causes a further decrease in global mean temperature of 0.34 °C. As in Stap et al. (2014), a shift in the mid-point of ocean circulation to 15° S produces surface temperatures which are in better agreement with LGM reconstructions. In particular, shifting the ocean circulation mid-point 10° further south decreases temperature and precipitation rates in the SH extratropics and increases temperature and precipitation in the NH extratropics, with negligible changes in global mean temperature and precipitation. Finally, reductions in the strength of overturning circulation cause simultaneous warming in the tropics and cooling in the higher latitudes of each hemisphere, which generate global mean cooling and a larger disagreement with the PMIP3–PMIP4 ensemble and LGM reconstruction from Annan et al. (2022).

Annan et al. (2022)

Table 4 Difference in zonal- and annual-mean surface air temperatures (°C) between the LGM and the PI. Inso refers to the standard ZEMBA-LGM experiment shown in Fig. 9.

Table 5 Difference in zonal- and annual-mean precipitation rates (mm d⁻¹) between the LGM and the PI. Inso refers to the standard ZEMBA-LGM experiment shown in Fig. 9.

4 Discussion

Our simulation of the pre-industrial period demonstrates that ZEMBA can describe the zonally averaged climate of this period with reasonable accuracy. Surface air temperatures are in strong agreement with NorESM2 and ERA5 1940–1970 both annually and seasonally (Fig. 2a–c), generally falling within 4 °C of those estimated by NorESM2 (Fig. 2d–f). Most notably, given the simplicity of ZEMBA, the model shows success in emulating precipitation and snowfall rates for the pre-industrial climate (Fig. 3a), particularly at the polar latitudes. Furthermore, the model captures the poleward enhancement of the surface and planetary albedo due to the presence of snow and sea ice cover (Fig. 4), with the seasonal range of sea ice cover in the NH corresponding nicely to both NorESM2 and ERA5, varying between ∼ 0.8 and ∼ 1.4 × 10¹³ m² (Fig. 5b). The ocean transport model shows success in emulating zonally averaged ocean heat transport, which makes a significant contribution to total heat transport at the lower latitudes, peaking with values of around 1 and 1.5 PW for the SH and NH, respectively (Fig. 7b). Similarly, the model predicts that atmospheric heat transport dominates in the middle to high latitudes and reaches its maximum levels at around ∼ 45° in each hemisphere. Finally, looking specifically at the partition of atmospheric heat transport into its dry static and latent components, the model captures the relative contributions of each flux to the total atmospheric heat transport, with dry static fluxes peaking at ∼ 15° N/S and latent heat transport peaking at 45° N/S.

Despite these favourable comparisons with state-of-the-art climate models and reanalysis data, there are limitations of ZEMBA in its current state. For example, while the model does reasonably well at capturing the radiative fluxes at the TOA and the surface (Fig. 6), one of the largest discrepancies between ZEMBA and NorESM2 or ERA5 1940–1970 resides in the outgoing longwave radiation at the TOA. As the outgoing longwave flux is overestimated around the tropics, the net radiation received at these latitudes is up to ∼ 25 W m⁻² lower than NorESM2 and ERA5 1940–1970 (Fig. 6c). Given the comparatively large surface area in the tropics, the underestimation of these net radiative fluxes results in a significant reduction in the surplus radiative energy which needs to be transported poleward. Consequently, both the total and the atmospheric heat transport in ZEMBA is much lower (Fig. 7a–b), with atmospheric heat transport peaking at ∼ 3 PW in each hemisphere, whereas it should be closer to 5 PW (Trenberth and Fasullo, 2017). In addition, one of the challenging aspects of a simple model such as ZEMBA is the underestimation of the winter maximum in snow cover over land (Fig. 5a), resulting in a smaller seasonal amplitude in snow cover. However, we note that the annual-mean rates of snowfall correspond well to NorESM2 and ERA5 1940–1970 over the northern high latitudes (Fig. 3a). Similarly to NH snow cover, the seasonal range of sea ice cover in the SH is reduced when compared to NorESM2 and ERA5 1940–1970 (Fig. 5b) and snowfall rates over Antarctica are underestimated by about 50 % (Fig. 3a).

Given their simplicity, the ability of EBMs like ZEMBA to accurately reproduce the latitudinal pattern of surface temperature is impressive but has been well established since the original works by Budyko (1969) and Sellers (1969). In reference to other studies, ZEMBA is most similar to later iterations of EBMs that contain vertical resolution (i.e. an atmospheric and surface layer), the division of the surface into land and ocean, and a seasonal cycle in insolation (Peng et al., 1987; Harvey, 1988; Bintanja, 1997). In comparison to the “present-day” simulation from the EBM developed by Bintanja (1997), which effectively serves as the basis for ZEMBA, most climate variables exhibit qualitative agreement with the pre-industrial output shown in Sect. 3.1, although the global mean surface air temperature reported in Bintanja (1997) is 1.09 °C warmer due to a higher 350 ppm CO₂ forcing in their model experiment. Other than the hydrological cycle, the primary contrast between the models lies in their representation of land-based snow cover. In the Bintanja (1997) study, the seasonal cycle in snow cover is parameterized as a function of surface air temperature and appears to be in better agreement with observations than ZEMBA, which explicitly calculates snow coverage over land. Overall, many such features of the Earth's climate are relatively well produced by various EBMs (Peng et al., 1987; Harvey, 1988; Bintanja, 1997), including ZEMBA, in part because they are “tuned” to match the present-day climate. Of greater importance for EBMs than achieving an exact replication of observations and/or GCM output is the investigation of different climate processes and feedbacks, which necessitates appropriate sensitivities to alterations in external or internal forcings.

The responses of ZEMBA to both a doubling of the atmospheric CO₂ concentration and a 2 % increase in solar insolation are in qualitative and quantitative agreement with one another (Fig. 8), including polar amplification of warming, which is stronger in the NH than the SH and concentrated in the winter months. The climate sensitivity of ZEMBA is larger than that of the EBM of Bintanja (1997) and similar such studies (Peng et al., 1987). For example, the S_o+2 % experiment generates much stronger annual-mean warming in the high northern latitudes (>14°) than Bintanja (1997) (<5°). In addition, the ECS of ZEMBA (3.6 °C) is higher than that of the original model of 1.9–2.2 °C (Bintanja, 1997; Stap et al., 2014). However, more recent EBMs have generated global mean and polar-amplified warming in response to CO₂ forcing which is similar to or even greater in magnitude than that of ZEMBA (Roe et al., 2015; Södergren et al., 2018; Feldl and Merlis, 2021). Moreover, the ECS of ZEMBA is consistent with the ECS range projected by GCMs that contributed to the Coupled Model Intercomparison Project Phase 5 (CMIP5) and CMIP6 (Flato et al., 2013; Zelinka et al., 2020) and other works (Collins et al., 2013). The strong seasonal asymmetry of surface warming in the polar regions (Fig. 8b–d) is in accordance with both observations (Screen and Simmonds, 2010) and GCM results (Holland and Bitz, 2003). As in Bintanja (1997), the winter maximum in surface air temperature is related to sea ice loss, causing a greater absorption of shortwave radiation, which is released from the ocean surface into the overlying atmosphere during the winter months. For the 2xCO₂ experiment, the polar amplification in warming – normalized (divided) by the global mean warming – exceeds 4 in the high latitudes of the Arctic and reaches up to 1.5 in the Antarctic (Fig. 8a). While polar amplification in warming that is strongest in the NH is in agreement with both observations and GCM results, warming in the high northern latitudes that reaches 4 times the global mean resides in the upper boundary of estimates from GCM simulations (Holland and Bitz, 2003; Hahn et al., 2021). When averaged from 60 to 90° N, however, the normalized Arctic warming is 3.2, which is still less than the 3.5 estimated over the historical period from 1979 to 2014 according to the observational dataset HadCRUT5 (Hahn et al., 2021). Overall, we see that ZEMBA exhibits a climate sensitivity – in terms of both global mean and polar-amplified warming – that appears broadly consistent with other EBMs, GCMs and observations.

There are comparatively few EBMs which incorporate a hydrological cycle (Jentsch, 1991; Kukla et al., 2023). Jentsch (1991) developed an EBM, consisting of a vertically averaged atmosphere overlying an ocean layer, to study the influence of the hydrological cycle on climate. Using an optimization procedure, the modelled precipitation and evaporation, amongst other climate variables, compare very well with contemporary observations (Jentsch, 1991). While both ZEMBA and the work of Jentsch (1991) contain a hydrological cycle, ZEMBA includes land cover and a seasonal cycle, which is suitable for investigating the response of climate to variations in the orbital parameters. Kukla et al. (2023) introduced a moist static energy balance model with a hydrological cycle coupled to a carbon cycle model, suited for studying the long-term relationship between the carbon cycle, hydrological cycle and climate. However, the modelled precipitation does not impact the surface albedo, which is instead simplified as a function of surface temperature, and, like the work of Jentsch (1991), it lacks the seasonal insolation cycle necessary to study the climate response to Milanković cycles. Previous EBMs used for studies of glacial–interglacial cycles prescribe either a present-day (Pollard, 1978) or a spatially uniform (Huybers and Tziperman, 2008) distribution of precipitation, which is not perturbed across climates, or instead parameterize precipitation and snowfall as a function of surface air temperatures and ice sheet size (Stap et al., 2014). The inclusion of a hydrological cycle enables precipitation and snowfall to be calculated internally in ZEMBA, which can account for the influence of changes in both local air temperatures and poleward moisture transport for precipitation and snowfall rates. It has been suggested that changes in the Earth's obliquity – by altering the meridional gradient in insolation – have a relatively strong influence on ice sheet volume due to changes in poleward moisture transport (Raymo and Nisancioglu, 2003; Nisancioglu, 2004). The inclusion of a hydrological cycle enables ZEMBA to explore these processes in the context of glacial–interglacial cycles.

It is important that ZEMBA, intended for studies of glacial–interglacial cycles, can simulate climates other than the present-day or pre-industrial periods. When our simulation of the LGM is compared to other reconstructions, the model compares favourably regarding changes in both surface air temperature (Fig. 9a) and precipitation (Fig. 9b). Indeed, ZEMBA captures the polar amplification of cooling in both hemispheres due to positive feedbacks relating to snow cover and sea ice expansion, although this cooling appears somewhat underestimated around Antarctica and slightly overestimated in the northern high latitudes. Moreover, the global mean cooling of −4.11 °C is similar to that estimated by the ensemble averages from PMIP3 and PMIP4 (Kageyama et al., 2021) and the data assimilation product from Annan et al. (2022). It should be noted that other data assimilation reconstructions suggest much stronger global cooling of between 6.1 and 6.8 °C (Tierney et al., 2020; Osman et al., 2021), which perhaps signals the importance of including feedbacks relating to clouds, dust and/or vegetation for reproducing the LGM cooling. However, the reconstructions from Tierney et al. (2020) and Osman et al. (2021) are based on a single climate model (CESM1.2), which produces one of the coldest LGM climates in PMIP4, whereas the Annan et al. (2022) reconstruction incorporates the wide range of climates generated across the PMIP ensemble.

Evaluating changes in meridional heat transport during the LGM is made difficult by the large (and often conflicting) range of total, atmospheric and ocean heat transport generated across PMIP3 and PMIP4 at all latitudes (Kageyama et al., 2021). For the ZEMBA-LGM experiment, changes in total heat transport are most intense in the NH (Fig. 9c), in response to large steepening of the Equator-to-pole temperature gradient. The upsurge in total heat transport between 15–45° N is somewhat consistent with PMIP3–PMIP4 (Kageyama et al., 2021), although most of this increase in ZEMBA is generated by stronger atmospheric heat transport (Fig. 9c–d), whereas PMIP3–PMIP4 also shows stronger ocean heat transport at these latitudes. Stronger ocean heat transport in PMIP3–PMIP4 can be attributed to a stronger and sometimes deeper Atlantic Meridional Overturning Circulation (AMOC) (Kageyama et al., 2021), which is inconsistent with proxy-based reconstructions of ocean circulation (Lynch-Stieglitz et al., 2007; Gebbie, 2014; Du et al., 2020) and is associated with stronger surface winds over the northern North Atlantic (Muglia and Schmittner, 2015). For ZEMBA, on the other hand, there is a large decreases in NH ocean heat transport during the LGM, reaching ca. −0.5 PW at 60° N (Fig. 9c). The prominent reduction in ocean heat transport, not observed in any of the PMIP3–PMIP4 simulations, is associated with the expansion of sea ice and a surface ocean layer that resides near the freezing point of seawater. Consequently, horizontal heat fluxes due to both advection (representing overturning) and diffusion (representing eddies and gyres) drop to zero at latitudes now covered by sea ice, which perhaps highlights the limitation of the simplified ocean model used in ZEMBA, where meridional heat transport is limited to the surface and bottom ocean layers.

The sensitivity experiments performed for the LGM boundary conditions (Fig. 10) recreate the experiment made by Bintanja and Oerlemans (1996) using the original EBM. For ZEMBA, the global mean cooling caused by the combined ice sheet and CO₂ forcing is 4.11 °C, which is stronger than the 3.3 °C of cooling noted by Bintanja and Oerlemans (1996) and is in better agreement with LGM reconstructions. The stronger cooling in ZEMBA can be partially attributed to differences in the LGM boundary conditions used in this study, which involves a stronger CO₂ reduction and larger LGM ice sheet extents than those in Bintanja and Oerlemans (1996). In our study, the addition of LGM ice sheets generates cooling and drying that are primarily localized to the NH, whereas CO₂ lowering generates widespread cooling and drying with stronger global mean cooling than that induced by ice sheet expansion (Tables 4 and 5). Additionally, by shifting the mid-point of ocean circulation from 5 to 15° S, surface air temperatures are in better agreement with the LGM reconstructions, as first noted by Stap et al. (2014). As for the strength of ocean circulation, Bintanja and Oerlemans (1996) found that reducing flow velocities in their ocean model led to temperatures in better agreement with contemporary LGM temperature reconstruction. In this study, however, reductions in ocean overturning instead lead to larger discrepancies with modern LGM constructions. Indeed, it remains uncertain to what degree overturning changed during glacial conditions, with some studies suggesting the mean state of overturning was not “sluggish” across glacial cycles (Bohm et al., 2015) or the LGM (Lynch-Stieglitz et al., 2007). Therefore, for simulations of glacial–interglacial cycles using ZEMBA, it may be important to consider the impact of variations in both the mid-point and the overall strength of ocean circulation.

ZEMBA is intended to be used as a computationally efficient tool for studies of the glacial–interglacial cycles of the Quaternary. The PI and the LGM experiments indicate ZEMBA is able to simulate glacial and interglacial climate states for a given insolation, CO₂ and ice sheet extent. As the model includes (1) both hemispheres, (2) a seasonal cycle and (3) a hydrological cycle, it is able to explore mechanisms invoked to explain the dearth of precession cycles in ice volume observed during the Early Pleistocene, such as (1) out-of-phase precession cycles between the hemispheres (Raymo et al., 2006), (2) a counterbalancing between summer insolation intensity and summer duration (Huybers, 2006), or (3) obliquity-induced variations in atmospheric moisture transport (Raymo and Nisancioglu, 2003). In future work, we intend to explore both the equilibrium of ZEMBA and its transient response to changes in the Earth's obliquity and precession, prior to simulations of the glacial–interglacial cycles of the Early Pleistocene via coupling to an ice sheet model.

5 Conclusions

In this study, a simple climate model (ZEMBA) is introduced to simulate zonally averaged climate fields including surface temperatures and precipitation. ZEMBA is largely built on the zonally averaged energy balance climate model from Bintanja (1997), comprising an atmospheric layer overlying a surface divided into a land component and a six-layered, zonally averaged ocean transport model. Unlike its predecessor, ZEMBA incorporates a hydrological cycle to estimate snowfall and precipitation with latitude.

Simulations of the pre-industrial period compare favourably with GCMs and reanalysis data, including surface temperatures (and their seasonal cycle), precipitation, surface and TOA radiative fluxes, sea ice, snow cover, and meridional heat transport. However, the underestimation of the net TOA radiation received in the tropics leads to a reduction in atmospheric heat transport, and there is also an underestimation of the seasonal amplitude in snow cover over land.

The responses of ZEMBA to increases in the atmospheric CO₂ concentration or the solar constant are in qualitative agreement with other EBMs, GCMs and observations, such as polar amplification in surface warming, which is strongest over the NH and focused in the winter months. The new additions to ZEMBA appear to increase climate sensitivity in comparison to older EBMs, but its results are still broadly consistent with the global mean and polar-amplified warming projected by climate models of higher complexity.

As the purpose of ZEMBA is for studies of the glacial–interglacial cycles of the Quaternary, it is important that the model can simulate climates other than present-day or pre-industrial climates. A simulation of the LGM indicates ZEMBA is able to capture changes in surface temperature and precipitation in qualitative and quantitative agreement with state-of-the-art climate models and data assimilation products, despite neglecting climate feedbacks relating to dust, vegetation and clouds. In particular, ZEMBA reproduces the polar amplification of cooling in both hemispheres and global mean cooling in accord with reconstructions from more elaborate models.

The overall conclusion from this study is that ZEMBA is suitable for studies of climatic change on large spatial and temporal scales, with a particular emphasis on glacial–interglacial cycles and the response of the climate system to changes in the orbital parameters. In future work, we intend to explore both the equilibrium of ZEMBA and its transient response to changes in the Earth's obliquity and precession, prior to simulations of the glacial–interglacial cycles of the Early Pleistocene via coupling to an ice sheet model.

Appendix A: Additional sensitivity experiments

A1 Sensitivity to prescribed cloud cover

To assess the limitations of prescribing a single cloud cover fraction from a pre-industrial simulation of NorESM2, we repeat the ZEMBA-PI simulation with different choices of zonal-mean cloud cover. In the ZEMBA-PI_CESM2 and ZEMBA-PI_MRI-ESM2 experiments, we force ZEMBA with pre-industrial cloud cover fractions taken from the Community Earth System Model 2 (CESM2) and the Meteorological Research Institute Earth System Model Version 2.0 (MRI-ESM2), respectively. CESM2 generates larger cloud cover fractions (ranging from 5 % to 25 %) relative to NorESM2 for the pre-industrial period, whereas cloud cover fractions from MRI-ESM2 and NorESM2 correspond more closely (Fig. A1a, c). In the ZEMBA-PI_ERA5 and ZEMBA-PI_CERES experiments, ZEMBA is forced with cloud cover fractions taken from the ERA5 atmospheric reanalysis averaged from 1940 to 1970 and from the Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled product averaged from 2005 to 2015, respectively. Similarly to CESM2, the CERES product contains much larger cloud cover fractions over the tropics and mid-latitudes but noticeably less cloud cover in the polar latitudes (Fig. A1a, c). On the other hand, the ERA5 dataset contains similar cloud cover to NorESM2 at the lower latitudes but has higher cloud cover in the polar regions (Fig. A1a, c).

Figure A1b shows the surface temperatures generated by ZEMBA when forced by these different cloud cover fractions, and Fig. A1d shows these anomalies relative to the standard PI simulation of ZEMBA (using NorESM2 PI cloud cover). When ZEMBA is forced by MRI-ESM2 and ERA5 cloud cover, which correspond closely to NorESM2 over most of the tropics and mid-latitudes, the differences in zonal-mean temperature are quite small. Changes in global mean temperature do not exceed 0.25 °C between these three PI simulations. However, when forced by CESM2 or CERES cloud cover fractions, which are generally much higher than NorESM2, there is a strong cooling effect. The decrease in temperature is unsurprising given that the net effect of more cloud cover is to cool the Earth for the radiation parameterization used in ZEMBA from Bintanja (1996). However, the global mean cooling for CERES (2.4 °C) is larger than that of CESM2 (1.67 °C), despite CESM2 having a consistently larger cloud cover. In particular, using CERES cloud cover generates much stronger cooling at high latitudes, in the region where CERES has very low cloud cover and where CESM2 has very high cloud cover. This suggests that, while the net global effect of clouds is to cool the Earth, the warming effect of clouds (via longwave radiation) outweighs the cooling effect (via shortwave radiation) in the polar regions. Indeed, ERA5 also generates warmer temperatures in the polar regions, where cloud cover fractions are higher than for NorESM2. Overall, we see that the choice of cloud cover fractions can have a strong impact on surface air temperature for the pre-industrial period. While the differences are small compared to MRI-ESM2 and ERA5 cloud cover fractions, they can become substantial for CESM2 and CERES. The choice of a different cloud cover fraction in ZEMBA would require a retuning of other model parameters to ensure the model simulates surface air temperatures with reasonable accuracy for the pre-industrial period.

Figure A1 (a) Different values of zonal-mean cloud cover including NorESM2 (ZEMBA-PI in grey), CESM2 (ZEMBA-PI_CESM2 in blue), MRI-ESM2 (ZEMBA-PI_MRI-ESM2 in green), CERES 2005–2015 (ZEMBA-PI_CERES in black) and ERA5 1940–1970 (ZEMBA-PI_ERA5 in red); (b) the zonal-mean surface air temperature simulated by ZEMBA in response to these different cloud cover fractions; (c) the differences in zonal-mean cloud cover relative to NorESM2 PI cloud cover; and (d) the differences in zonal-mean surface air temperature (caused by using different cloud cover fractions) relative to the standard ZEMBA-PI simulation. Pre-industrial cloud cover from CERES and MRI-ESM2 is taken from the Earth System Grid Federation at https://esgf-node.llnl.gov/search/cmip6/ (last access: 25 October 2024). ERA5 cloud cover is taken from Hersbach et al. (2023) (https://doi.org/10.24381/cds.f17050d7, last access: 25 October 2024), and CERES is taken from https://ceres.larc.nasa.gov/data/ (last accessed: 25 October 2024).

A2 Sensitivity to key model parameters

We examine the sensitivity of ZEMBA to internal model parameters by replicating selected experiments from Bintanja (1997). In addition to the CO₂ (ZEMBA-2×CO₂) and solar-constant (ZEMBA-S_o+2 %) experiments described earlier, we perform simulations with perturbations in cloud cover amount (+20 %), the cloud optical depth (τ+1.4), the diffusion coefficients for both horizontal (D_o×2) and vertical (D_z×2) ocean heat transport, the turbulent heat flux coefficient (κ×2), and sea ice thickness (d_si×2). We also perform a new experiment with a perturbation in a Hadley cell constant (λ+0.6). Figure A2 shows the changes in global mean surface air temperature, global mean ocean temperature and the average Equator-to-pole temperature gradient driven by these large perturbations in model parameters.

The response of ZEMBA to perturbations in these model parameters is qualitatively similar to the responses reported by Bintanja (1997), though ZEMBA shows heightened sensitivity, as identified previously in the ZEMBA-2×CO₂ and ZEMBA-S_o+2 % experiments. This increased climate sensitivity in ZEMBA may be attributed to changes such as the parameterization of land surface albedo and the use of atmospheric heat transport set proportional to gradients in moist static energy, rather than temperature, which can enhance polar-amplified warming. Both the 2 % increase in the solar constant and the CO₂ doubling yield similar outcomes: global mean and polar-amplified warming leading to a reduced Equator-to-pole temperature gradient.

Perturbations in the cloud cover parameters cause the most significant changes in surface temperature. The net effect of increased cloud cover is to drive global mean cooling. Moreover, an increase in cloud optical depth (τ), which increases the cloud albedo, produces an even stronger cooling effect than changes in total cloud amount. As noted by Bintanja (1997), the stronger sensitivity to an increased cloud optical depth is due to its sole effect of enhancing the reflection of incoming shortwave radiation at the TOA. In contrast, when cloud cover is increased, the enhanced shortwave reflection is partially offset by decreased outgoing longwave radiation, which moderates the overall cooling effect.

Adjustments to the ocean heat transport coefficients, turbulent heat flux coefficients or sea ice thickness have more modest impacts on global mean temperature compared to cloud cover. Doubling the eddy and gyre diffusion coefficient (D_o) marginally raises polar temperatures without significantly affecting the global mean. Similarly, doubling the vertical diffusion coefficient (D_z) has a negligible effect on air temperature but markedly affects mean ocean temperature. It should be noted, however, that a large proportion of ocean heat transport is carried out by the prescribed ocean overturning, which is unaffected by these diffusion coefficients. Increasing κ enhances the surface-to-atmosphere heat fluxes, which leads to surface cooling and thereby sea ice and snow expansion and ultimately results in global mean cooling. Increasing d_si results in a slight global temperature decrease by reducing the seasonal variability in sea ice, thereby increasing its extent during summer months (Bintanja, 1997).

Figure A2 Anomalies in global mean surface air temperature, global mean ocean temperature, and the air temperature difference between the equatorial (0–10°) and polar regions (80–90°) for changes in the solar constant (S_o), atmospheric CO₂ level (CO₂), cloud amount, cloud optical depth (τ), ocean diffusion coefficient for horizontal eddy and gyre heat transport (D_o), ocean diffusion coefficient for vertical heat transport (D_z), coefficient for turbulent heat fluxes (κ), sea ice thickness (d_si) and a Hadley cell parameter (λ).

Finally, to evaluate sensitivity related to the Hadley cell parameterization, we modify λ – representing the fractional difference between the upper branch's uniform moist static energy and the surface moist static energy at the Equator (${\stackrel{\mathrm{‾}}{m}}_{\mathrm{eq}}$). In the original formulation from Siler et al. (2018), λ=1.06, which we adjusted to 1.03 for an improved simulation of PI precipitation. In this sensitivity experiment, we increase λ to 1.09. While changing λ affects zonal precipitation in the tropics, the total atmospheric energy transport remains governed by meridional gradients in moist static energy, leaving global mean temperature largely unaffected. Additionally, assumptions involving the Ω weighting function, which determines the Hadley cell's dominance in heat transport, were explored by Siler et al. (2018), who noted that varying representations of Ω did not alter the primary climate response.

Appendix B: List of model parameters

In Table B1 we provide definitions for all the model parameters listed in Sect. 2.

Table B1 All model parameters listed in the atmospheric, land and ocean components of the model.