TransEBM v. 1.0: Description, tuning, and validation of a transient model of the Earth’s energy balance in two dimensions

. Modeling the long-term transient evolution of climate remains a technical and scientiﬁc challenge. However, understanding and improved modeling of the long-term behavior of the climate system increases conﬁdence in projected changes in the mid-to long-term future. Energy balance models (EBMs) provide simpliﬁed and computationally efﬁcient descriptions of long timescales and allow large ensemble runs by parameterizing energy ﬂuxes. This way, they can be used to pinpoint periods 5 and phenomena of interest. Here, we present TransEBM, an extended version of the two-dimensional energy balance model by Zhuang et al. (2017a). Transient CO 2 , solar insolation, orbital conﬁguration, ﬁxed ice coverage and land-sea distribution are implemented as effective radiative forcings at the land surface. We show that the model is most sensitive to changes in CO 2 and ice distribution, but the obliquity and land-sea mask have signiﬁcant inﬂuence on modeled temperatures as well. We tune TransEBM to reproduce the 1960 – 1989 C.E. global mean temperature and the equator-to-pole and seasonal temperature 10 gradients of ERA 20 CM reanalysis (Hersbach et al., 2015). The resulting latitudinal and seasonal temperature distributions agree well with reanalysis and the general circulation model (GCM) HadCM 3 for a simulation of the past millennium (Bühler et al., 2020). TransEBM does not represent the internal variability of the ocean-atmosphere system, but non-deterministic elements and non-linearity can be introduced through model restarts and randomized forcing. As the model facilitates long transient simulations, we envisage its use in exploratory studies of stochastic forcing and perturbed parameterizations, thus 15 complementing studies with comprehensive GCMs


Introduction
Observational records of climate and proxies are limited in both time and space (Schmidt et al., 2014;Hersbach et al., 2015). As such, climate models are vital for climate research since they can fill gaps left by proxy and observational records, for example with respect to transitions and hysteresis in climate (e.g., Rahmstorf et al., 2005;Abe-Ouchi et al., 2013;Ghil, 2015;Fraedrich et al., 2016;Stap et al., 2017), or variability in the climate system (Huybers and Curry, 2006;Laepple and Huybers, 2014; that an increase in variability indicates a predisposition of Earth's climate towards tipping into a different state in response to small changes in boundary conditions. Ghil (2015) investigated climate variability and sensitivity from a system dynamics standpoint to understand variability on all timescales. For this, the study used zero-to two-dimensional EBMs as part of a model hierarchy to show how multiple equilibria arise and how internal variability interacts with forcings. Dortmans et al. (2019) also found evidence for bifurcations and possible transitions between stable states in relation to glaciation. They connected these 70 transitions to particular points in the geological past using a zero-dimensional EBM. At these points, they investigated possible mechanisms that explain sudden glaciations or unusually warm periods, focusing mainly on CO 2 , water vapor, as well as oceanic and atmospheric heat transport.
Further examples of studies of Earth's climatic history include the simulation of a possible climate on the supercontinent Pangaea during the early Late Permian , or on Gondwana (Crowley et al., 1986(Crowley et al., , 1987. Hyde et al. (1989) 75 simulated the past 18, 000 years with snapshots taken every 3000 years and Short et al. (1991) modeled the past 800, 000 years in selected regions using snapshot simulations in response to orbital forcing. For such long-term studies in particular, transient simulations provide insights into the build-up of forcings, feedbacks, and their interactions in the climate system. With respect to transient simulations, Harvey and Schneider (1985) presented a pioneering study, which used a hemispherically averaged box-type EBM to examine effects of solar, CO 2 , and volcanic forcing. More recently, Stap et al. (2017) created transient 80 simulations of the past 38 million years to study the impact of ice sheets on the relationship between temperature and CO 2 by coupling an one-dimensional EBM to an one-dimensional ice sheet model. Moreover, EBMs have been applied to study scaling of the power spectrum of temperature, for example by Rypdal et al. (2015), as well as the temperature response to stochastic forcing (e.g., Kim and North, 1991;Rypdal et al., 2015). Fraedrich et al. (2016) and Meyer et al. (2018) used them to analyze memory effects in the climate system. Dommenget et al. (2019) 85 built a large database of more than 1, 300 experiments with a two-dimensional EBM. They explored how different feedbacks and components in the climate system such as the ice-albedo feedback, clouds, or the hydrological cycle respond to changes in forcing. Many of these studies provide robust results only when using large ensembles of runs that represent many possible realizations of climate or explore the parameter space of the models.
We present TransEBM v. 1.0, a transient extension of the model by Zhuang et al. (2017a), a two-dimensional EBM that 90 simulates the surface temperature field based on the Earth's energy balance. Briefly, the model works as follows: Local albedo determines how much of the incoming radiation is reflected back into space by both atmosphere and surface. The atmosphere is averaged with no local differences and only the surface temperature is computed. The surface temperature changes according to the absorbed energy and the surface's heat capacity. The energy that enters the system is transported so as to balance the Figure 1. Land, sea, and ice configuration with T42 resolution for 1950 CE as used in modern model runs and given by Zhuang et al. (2017a).
Four surface types are differentiated: ocean, land surface, sea-ice, and land-based ice/ snow. temperature gradients caused by latitudinal differences in incoming radiation and the inhomogeneous distribution of land and 95 ice masses. Earth then emits longwave radiation depending on the local surface temperature.
The model is presented in detail in Sec. 2, which covers the model description (Sec. 2.1), the parameterization and forcings (Sec. 2.2), the code structure and benchmarking (Sec. 2.3), and the parameter tuning (Sec. 2.4). Section 3 validates the model by starting with the evaluation of the tuning using reanalysis data (Sec. 3.1). Furthermore, it presents a sensitivity test with respect to the forcings (Sec. 3.2) and the results of a past millennium simulation, which are compared to GCM simulations and 100 a multi-proxy reconstruction (Sec. 3.3). We present a short tutorial in Sec. 4, followed by a discussion of the model's strengths and weaknesses, as well as potential applications and extensions in Sec. 5.

Model description
TransEBM computes a surface temperature field T (r, t) that depends on time t and two spatial dimensions in latitude and 105 longituder = (θ, φ) according to with heat capacity C(r), outgoing long-wave radiation parameters A and B, heat transport coefficient D(r), solar constant S 0 , insolation forcing S F (r, t), and albedo a(r). As detailed in Zhuang et al. (2017a), this elliptic, partial differential equation is solved for 48 time steps per year using the full multigrid method, which dampens both high-and low-frequency errors (Briggs Figure 2. Conceptual view of the forcings, energy transport, and parameterizations in TransEBM. Solar insolation S0, the orbital parameters, -eccentricity, obliquity, precession of the perihelion -CO2, land, sea, and ice distribution serve as forcings. The surface types differ by albedo a, which affects the insolation field, and heat capacity C. The energy from incoming solar radiation is transported poleward according to the thermal conductivities κ. The outgoing radiation is calculated based on the parameters A, B, and the response to CO2 scaled by CO2,s.
snow -as shown in Fig. 1

Forcings and parameterization
Four broad categories of parameters exist: heat capacity, albedo, thermal conductivity, and outgoing radiation. Figure 2 shows a conceptualized view of the energy transport of the model with all model parameters and forcings.
The simulated heat capacity changes with the surface type as listed in Table 1. Zhuang et al. (2017a) base the oceanic heat capacity on the assumption of a mixed layer of about 70 m. This can be changed by tuning the oceanic heat capacity. The 125 tuned heat capacity of TransEBM corresponds to a mixed layer of about 96 m. The contribution of the deep ocean that impacts timescales of order 10 2 years and higher (Rohling et al., 2012) is not modeled. It is possible to differentiate further surface types than the ones currently used in the model. In fact, additional surface types to distinguish between water reservoirs like the Pacific, Atlantic, and Indian ocean, as well as the Mediterranean and inland seas exist for the parameterization of heat capacity.
Using them, the impact of the deep ocean could be modeled more accurately. However, other parts of the parameterization like 130 Similarly, the albedo depends on the geography, gets initialized in the beginning of a model run, and then remains constant. latitude-dependent. Figure 3 shows the albedo distribution of a present-day continental configuration as seen in Fig. 1. The highest albedo values are in snow-covered areas and the lowest close to the equator.  Heat transport due to sensible and latent heat, as well as oceanic heat transport across isotherms is described using a heat transport coefficient (κ land − κ land,NP ) sin 5 (θ) + κ land,NP , Northern Hemisphere (κ land − κ land,SP ) sin 5 (θ) + κ land,SP , Southern Hemisphere (2) 140 which is constant in time. Figure 2 presents a conceptual view of the heat transport in the model. Oceanic heat transport is distinguished from heat transport over land and ice-covered areas. Over land and ice-covered areas, the Northern and Southern Hemisphere are treated separately so as to reproduce differences in observed temperature patterns due to the asymmetric distribution and sizes of the land masses. This difference increases towards higher latitudes. Figure 4 shows the resulting latitude-dependent pattern of the heat transport coefficient. Heat transport is large close to the equator and declines towards 145 higher latitudes, where it evens out.
The parameterization of outgoing radiation as A + B · T follows Budyko (1968) and uses phenomenological coefficients A and B. Both, A and B, were estimated from satellite observations in previous studies (cf. Kyle et al., 1995;North and Kim, 2017). A captures the direct radiative effect of atmospheric CO 2 on radiative forcing. Following Arrhenius (1896), Myhre et al. Oceanic heat transport is treated separately from that over land and ice-covered areas. For the latter, Northern and Southern Hemisphere are separated with larger heat transport coefficients in the north.   (Kyle et al., 1995;155 North and Kim, 2017). Large uncertainties are attached to the measurements of both A and B (North and Kim, 2017). Table 4 summarizes the parameterization of outgoing radiation.
Insolation is the only component of the parameterization that changes for each time step with the seasonal cycle. It depends on the latitude-and time-dependent insolation S(θ, t) as well as the albedo a(r) via S F = S(θ, t)(1 − a(r)).
Here, the orbital configuration is given by declination δ, absolute hour angle of the Sun at sunrise and sunset H 0 , and the normalized Earth-Sun distance ρ which depends on eccentricity e and the positional angle of Earth on its orbit around the Sun ν as given counterclockwise from perihelion. Figure 5 shows the relationship between orbital forcing, insolation, and the modeled temperature response. For the three orbital forcings -eccentricity, obliquity, and precession -it shows the difference in insolation and temperature between the 170 possible extremes of the orbital forcings as listed in Table 5 during northern summer. When examining one orbital parameter, the other ones were set to their listed average values. The simulations were run in equilibrium mode with the tuned model (cf. Sec. 2.4) and the land-sea-mask as shown in Fig. 1, constant CO 2 = 315 ppm and S 0 = 1362 W m −2 .
Precession of the perihelion has the largest effect on the distribution of solar insolation and the resulting changes in temperatures, followed by obliquity. The effect of the changing eccentricity is the smallest. Larger eccentricity results in a global 175 decrease in temperatures between May and July when perihelion is set to occur at vernal equinox in March. This means that the Earth will be farther from the Sun in northern summer with higher eccentricities. The effect is larger in the Northern Hemisphere, particularly in mid-latitudes.
Between May and July, a large value of obliquity decreases insolation in the Southern Hemisphere and leads to an increase in the Northern Hemisphere. In the Southern Hemisphere, this effect is strongest in the mid-latitudes, but disappears over 180 high latitudes. In the Northern Hemispheres, insolation increases up to around 60 • N, and then is slightly smaller in the high latitudes, leading to the the largest increase in temperatures at high latitudes and over land in mid-latitudes. For all orbital parameters, the insolation and temperature patterns even out over the ice-covered regions near the poles.
For the high precession simulation, summer solstice and perihelion coincide, leading to increased insolation and temperatures between May and July. Overall, the strongest increase in insolation occurs in low northern latitudes and the highest increase in 185 temperatures over Asia. Temperatures grow especially large over land and in the northern high latitudes. Appendix B examines the effect of the orbital parameters on insolation and temperatures further by repeating this analysis for November to January and annual averages. Table 5. Orbital parameter extremes used to examine the maximal influence that the parameters have on the modeled distributions of insolation and temperature as shown in Fig. 5. To examine the effect of one parameter, the other two were set to average values and an equilibrium simulation with land, sea, and ice configuration as in Fig. 1

Model extension and benchmarking
Based on the model code by Zhuang et al. (2017b), the model formulation was extended to enable transient and parallel runs, 190 as well as transient, user-defined forcings for CO 2 , solar insolation and orbital configuration. The land, sea, and ice configuration can be changed within a simulation via model restarts at prescribed intervals. Additionally, the model configuration was overhauled so as to allow as much customization for experiments as possible. Table 6 summarizes the major expansions of TransEBM v. 1.0 with respect to the model by Zhuang et al. (2017a). A default configuration is provided with the model code.
It contains detailed explanations for each variable that expand on the short tutorial provided in Sec. 4.

195
The model is written in Fortran 90 and uses bash scripts and CMAKE and as such should run on standard Linux or Unix machines. It contains a preprocessing component visualized in Fig. 6 that imports the experiment's configuration, prepares the geography, and computes the albedo field. Based on this, the main component sets up the specified experiment, key steps of which are presented in Fig. 7. The initial setup is then integrated until either equilibrium or the user-defined number of years is reached. Finally, the output is written to netCDF files.

200
Only the land, sea, and ice distribution -described in the map section of the configuration -is constant throughout a model run, all other forcings can be transient. For long-term simulations, changing it is essential. This can be achieved by restarting Figure 5. Difference that high and low orbital parameters as listed in Table 5 cause in solar insolation (left) and temperature field (middle) for the averages of the months May, June, and July (MJJ). A schematic illustration of the orbital parameters is added to the right. Equilibrium simulations with land, sea, and ice configuration as in Fig

Tuning methodology
The model is tuned to the 1960 − 89 climatology from ERA20CM reanalysis (Hersbach et al., 2015). Target  To get a complete picture of how the parameters affect the modeled temperatures, parameter values that disagree with observations or physics were explored as well. During the actual tuning, however, only configurations that agree with observations and physics were considered.
Albedo changes have a strong effect on temperature means and extremes and as such affect all metrics as seen in Fig. 8.

220
They do not influence the shapes of the seasonal temperature distribution, though, since the albedo field is constant throughout the model year and as such affects the insolation at all time steps the same way. The thermal conductivities κ describing the heat transport alter primarily the seasonality and equator-to-pole gradient. Decreasing them sharpens the latitudinal distribution and increases the seasonality and vice versa for increased thermal conductivities. In comparison to albedo, they have a small 12 https://doi.org/10.5194/gmd-2020-237 Preprint. Discussion started: 9 October 2020 c Author(s) 2020. CC BY 4.0 License. effect on global mean temperature and produce differences of about 5 • C over the whole explored parameter space. Changing 225 the heat capacities affects the equilibration time and the seasonal cycle, both in regard to the lag of the temperature response to the forcing and the absolute temperature difference between warmest and coldest month. Figure 9 shows that the outgoing 13 https://doi.org/10.5194/gmd-2020-237 Preprint. Discussion started: 9 October 2020 c Author(s) 2020. CC BY 4.0 License.  out affecting the other metrics. Minimizing the root mean square errors was prioritized over agreement with the GM T from ERA20CM GM T ERA = 13.27 • C. With respect to GM T , any tuning options were only discarded if they produced a change in GM T by several degrees. Since RM SE lat was the highest, with significant disagreement observed in the polar regions, it was the initial focus of the tuning. Table 8 shows the metrics before and after tuning. The agreement between the latitudinal dis- Both A and B affect all tuning metrics, whereas changes in CO2,s barely make a difference.
tributions in ERA20CM and TransEBM and the seasonal distributions in the Southern Hemisphere improved. RM SE seas,nh 235 increased slightly, but since this produced a large improvement in the other two RM SEs, this increase was considered acceptable. Although GM T was barely considered during the tuning, the final GM T is closer to the target. Table 7. List of the parameter groups, namely heat capacities C, albedos a, thermal conductivities κ, and outgoing radiation parameters A ref , B, and CO2,s, and which of the tuning metrics GM T , RM SE lat , and RM SEseas they affect based on the results shown in Fig. 8 and 9.
affects parameter GM T RM SE lat RM SEseas The model was tuned using the 1960 -1989 climatology of a simulation of the twentieth century. As a first validation of that 240 tuning, we compare that whole simulation to ERA20CM reanalysis data (Hersbach et al., 2015) and a simulation with the EBM based on the parameterization by Zhuang et al. (2017a). Both simulations used the constant modern land, sea, and ice configuration shown in Fig. 1, with CO 2 data from Köhler et al. (2017) and S 0 values from Coddington et al. (2016). mechanisms that enhance the temperature response to increasing CO 2 , leading to the EBMs underestimating greenhouse gas related warming. Also unlike the reanalysis, they show smoother changes and less inter-annual variability. The eleven-year solar cycle, which was superimposed, dominates inter-annual changes in GMT in the simulations. fitted to the seasonal cycle for every grid box. The phase lag is then the difference between the occurrence of the extrema in the fit to the extreme in the insolation on the 21st of June and December. Figure 11 shows the resulting phase lags for TransEBM and ERA20CM. Areas where the seasonal cycle does not follow a sine were excluded from the analysis. Appendix C provides details on the exclusion process. The resulting pattern in phase lags highlights that TransEBM produces larger lags. Over the oceans, the modeled phase lag is significantly higher than that in the reanalysis. Over land and ice, TransEBM simulates longer 260 phase lags as well, but the difference to the reanalysis is smaller than in oceanic areas.

Sensitivity test
To test TransEBM's response to different forcings, we designed a three-step sensitivity test. The model is initialized for condi-270 tions of the Last Glacial Maximum (LGM) with the forcings listed in Table 9. This simulation is run for 100 years with constant forcings. Then, either one forcing or all of them shift to those found during the mid-Holocene warm period around 6000 before present (BP). After another 100 years, the simulation switches to conditions for the year 2000 CE. Overall, the simulations thus cover three states in Earth's history with distinct climatic conditions that span a large range of the natural variability in forcings.
Only for the solar constant there would be no significant change, which is why artificial steps were introduced that reveal its  Fig. 1 and App. A1a, respectively. The map for the mid-Holocene warm period was then the result of an interpolation between the two according to the sea level reconstruction by Grant et al. (2012) using the bathymetry by the National Geophysical Data Center (1993) (cf. Fig. A1b).
280 Figure 13 shows the temperature time series resulting from the sensitivity test. Among the forcings, the distribution of ice sheets and sea ice and CO 2 have the largest effect on GMT. Despite its artificial enhancement, the solar constant has a smaller influence. Changes to the land-sea mask and to obliquity have a significant effect as well, whereas the variations due to eccentricity and precession are negligible. The sum of the temperature changes due to the individual forcing steps is slightly larger than the simultaneous change in all forcings, suggesting that the model contains some non-linear behavior.

Past millennium climate
To further examine TransEBM's abilities, we simulate the climate of the past millennium. The solar forcing includes volcanic eruptions and follows Neukom et al. (2019), as does CO 2 . The orbital parameters are calculated according to Berger (1978) and the land, sea, and ice distribution is that from Fig. 1. To ensure adequate time for spin-up, the simulation is initialised in 800 CE, but evaluated from 850 to 2000 CE. Figure 14   2. For using a compiler other than gfortan-7, change the compiler name in files config/config_module.sh, preprocess/preprocess.sh, and src/Makefile. as FC=gfortran-7. Additionally, the configuration folder needs to be provided as I/home/user/TransEBM/config.

Set the orbital parameters
For constant values, set ECC, OB, and PER. Or have them be calculated following Berger (1978). In the latter case, 345 provide a start year of the experiment as INITIAL_YEAR and set ORB_SRC="code" and ORB_FORCED=TRUE.

Define the length of the experiment
Provide the initial year as a common era value in INITIAL_YEAR, for example as 1950 or −21000. Set the number of years for which the simulation should be run as RUN_YEARS.

Choose, whether it is an equilibrium or transient simulation
For USE_EQUI=true, the experiment will be stopped once equilibrium or the number of RUN_YEARS is reached, whichever happens first. For USE_EQUI=false, the simulation will run for as many years as provided in RUN_YEARS.
7. Set the working and the output directory.
Set them using WRK_DIR and OUT_DIR.

355
Finally, the experiment should receive an ID.
Other possible changes to experiments include the option to change the parameterization and the writing of additional output.
Once the configuration is set up, the model can be run by executing run_ebm.sh with the configuration file as input: bash run_ebm.sh config/test.conf If no configuration is specified, the code defaults to config/default_config.conf.

Discussion
Building on previous modeling efforts from Budyko (1968) to Zhuang et al. (2017a), we have presented here TransEBM v. 1.0, a transient energy balance model that simulates Earth's temperature field based on CO 2 greenhouse forcing, orbital parameters and solar constant, as well as land, sea, and ice distribution. The provided tuning is in reference to the 1960-89 climatology from the ERA20CM reanalysis (Hersbach et al., 2015). This tuning improves the simulated seasonal cycle in the Southern

365
Hemisphere and especially the latitudinal temperature profile. Thereby, it improves the fit of the temperature profile in polar regions, showing that TransEBM can be used to model temperatures in ice-covered regions and as such provide ensemble runs exploring long timescales.
There remain small disagreements with reanalysis and GCM simulations over Antarctica, where TransEBM underestimates temperature at the border of the ice sheet and overestimates them towards the pole. This is due to the lack of topography in 370 the model, which is why the elevation of ice sheets cannot be considered. In the latitudinal temperature distribution, the dip in temperatures near the equator related to the Intertropical Convergence Zone (ITCZ) is not reproduced. This band of clouds and its influence on the atmospheric conditions around the equator cannot be parameterized by TransEBM with the latitudinal difference in albedo and heat transport parameters that remains constant across latitudes. Since atmospheric dynamics are reduced to poleward transport of heat across latitudes with no ocean currents or precipitation considered, the ITCZ is missing 375 from simulations. Therefore, TransEBM provides an incomplete picture of changes in the tropics, which strongly depend on this component of the hydrological cycle.
TransEBM follows the overall trend in GMT, but is unable to capture inter-annual variability due to the lack of dynamical processes modeled. With respect to the seasonal cycle, TransEBM lags the reanalysis data and produces smaller seasonal amplitudes. Moreover, there is less variability found in the spatial distribution of the amplitudes in the model. This mirrors the 380 findings of Hyde et al. (1989). The spatial structures of the amplitude and phase lag of the seasonal cycle that they found in their EBM are similar to the ones presented here (cf. Fig. 12 and 11). They also showed reduced spatial variability in comparison to observations. Hyde et al. (1989) found similar magnitudes of the phase lag (cf. Fig. 11 here and Fig. 3a in Hyde et al. (1989)), but the modeled amplitude differs (cf. Fig. 12 here and Fig. 1a in Hyde et al. (1989)). TransEBM produces amplitudes roughly twice as high as Hyde et al. (1989), in better agreement with the reanalysis. The seasonal cycle in TransEBM is a product 385 only of the seasonality in its insolation fields. As such, it lacks seasonality in the sea ice response, which would influence the parameterization of albedo, heat transport, and heat capacity as simulated by the model. First trials of including such a seasonal sea ice cycle that could lead to the inclusion of dynamical ice sheets have been unsuccessful since they destabilize the numerical scheme. Such a change to the model, therefore requires further research. Furthermore, the model does currently not explicitly consider energy conservation. This would be especially relevant for interactive changes to the land, sea, and ice 390 distribution between model restarts.
Overall, TransEBM makes large ensemble runs and long timescales accessible for exploration and can be used to design experiments with models more constrained by computational costs since TransEBM runs quickly with limited computational resources. As an easy-to-use climate model, it is suitable for use in education, and could be used to reproduce past studies with energy balance models for which no documented code exists (e.g., Hyde et al., 1989). In order to reproduce studies of hysteresis 395 and abrupt transitions between ice-covered and ice-free states (cf. North et al. (1983); Mengel et al. (1988); Crowley et al. (1989); Short et al. (1991)), an extension of the model with interactive ice sheets and sea ice would be necessary. When it comes to forcings, it underlines the sensitivity to ice configuration and CO 2 forcing of simulated climates in comparison to other forcings, which have smaller impacts on temperatures. Furthermore, it illustrates the large impact that the parameterization of the outgoing radiation as well as the albedo fields of ocean and land masses have on the model results.

Conclusions
We describe a transient energy balance framework suitable for simulating planetary temperatures on long timescales due to its computational efficiency. It can capture trends in the global mean temperature and simulates latitudinal temperature profiles that agree well with ERA20CM reanalysis data (Hersbach et al., 2015) and HadCM3 (Gordon et al., 2000;Pope et al., 2000), a general circulation model. There are differences with respect to the seasonal cycle especially in the Northern Hemisphere 405 related to the tuning of the model provided here. TransEBM v. 1.0 has a larger phase lag than that found in reanalysis and the amplitude of the seasonal cycle in the Northern Hemisphere is smaller. Overall, it simulates less variable temperature fields than either reanalysis, multi-proxy reconstructions, or more complex models due to its reduced complexity.
Among the modeled forcings, changes to the ice sheets, CO 2 level, and the sea ice configuration are shown to influence the simulated climate the most. Among the parameters, the albedo fields of the oceans and land masses as well as the parameter-410 ization of outgoing radiation have the largest impact and even small changes to these parameters can shift the model results significantly. With respect to the seasonality, the heat capacity of the oceans also carries some weight. The model behaves largely linearly, with only a small difference between the separate and simultaneous application of the forcings.
Desirable future improvements include the addition of a seasonal cycle in the sea ice, which currently remains constant during a model run. The inclusion of this seasonality is likely to improve the modeled seasonal cycle. Such a change to the 415 model could encompass all ice masses and have them respond dynamically to the modeled temperature changes. For this to be physically reasonable, consideration of energy conservation within the model is necessary, which is especially relevant for the simulation of climate on long timescales with significant changes to Earth's ice masses. Overall, TransEBM v. 1.0 is accessible for large and long ensemble runs and thus can complement complex models, in particular to study the influence of orbital cycles as well as other radiative forcings on Earth's energy balance.

420
Code and data availability. The current version of the model is available at https://github.com/paleovar/TransEBM and licensed under the GNU General Public License Version 3. The version of the model used for the analyses presented here is archived on Zenodo (Ziegler and Rehfeld, 2020), together with all input data, configuration files, and scripts to run the model. The model simulations and all other data to produce the plots of this paper as well as the code that generated the plots is archived there as well.
Appendix A: Land, sea, and ice distributions used for the sensitivity test 425 Figure A1 shows the land, sea, and ice configuration used for the sensitivity test presented in Sec. 3.2. Panel (a) shows the map for the LGM taken from Zhuang et al. (2017a). Panel (b) shows a map corresponding to the state 6000 years BP. For this, the sea levels from Grant et al. (2012) were taken and applied to the bathymetry by the National Geophysical Data Center (1993) to obtain a land-sea distribution. Based on the difference in sea levels between this time and both the LGM and modern configurations from Fig. A1a and 1, respectively, the difference in ice cover was interpolated. The amount of ice corresponding 430 to the difference in sea level was removed from the LGM map according to a weighting. The weighting preferred the removal of ice closer to the equator and at the borders of ice cover, that is ice with more ice-free neighbours. Ice was removed randomly if only a fraction of the ice with the same weighting needed to be removed. This interpolation is flawed -it ignores the height of ice sheets, energy conservation, and many more aspects -but for the idealized sensitivity test it suffices.
Appendix B: Effect of the extremes in orbital forcing on November-December-January and yearly insolation and  All three orbital parameters show a pattern mirrored between the hemispheres from that between May and July shown in Fig. 5. For eccentricity this means that a high eccentricity corresponds to larger insolation and temperatures than are found for a lower eccentricity if obliquity and perihelion are set to the averages noted in Table 5. This net positive effect is larger in 445 the Southern Hemisphere with insolation growing towards middle southern latitudes. The largest increase in temperatures are found over Antarctica, followed by eastern Africa and the Middle East. Altogether, temperatures increase the least over the oceans.
A high obliquity causes a decrease in insolation in the Northern and an increase in the Southern Hemisphere in comparison to the low obliquity simulation. The largest decrease can be found in the middle northern latitudes, and the largest increase 450 occurs around 60 • S. With regards to temperatures, the largest increase occurs in the southern polar regions, but increased temperatures are simulated for the northern polar regions as well. The largest decreases are found over Asia and near the equator.
With the perihelion at the June solstice, insolation is lower between November and January in comparison to the perihelion occurring at the December solstice. This decrease is the largest in lower latitudes of the Southern Hemisphere. It translates to 455 decreased temperatures over land masses, Asia and Africa in particular, and towards the polar regions, which are affected by the lack of energy entering near the equator, which, if present, would be transported polewards. Precession acts strongest in lower southern latitudes, whereas obliquity and eccentricity affect the middle and high latitudes more. Figure B1. Difference in solar insolation (left) and temperature (right) caused by extremes in the orbital parameters as listed in Table 5  Repeating the analysis for yearly averages yields the results shown in Fig. B2. Most of the differences average out between the simulations with orbital extremes. For the precession of the perihelion in particular, there is virtually no difference between 460 the experiments, as panels (e) and (f) show. Both eccentricity and obliquity show a symmetric pattern in the annual averages.
The difference between obliquity extremes is the largest, confirming the results from the sensitivity study from Sec. 3.2, which also found that obliquity had the most notable effect on the global mean temperature among the orbital parameters. Obliquity affects the high latitudes the most, where a high obliquity increases insolation and temperatures. Close to the equator, there is a net negative effect on the insolation and temperatures. A high eccentricity, on the other hand, produces the largest effect close 465 to the equator, with both insolation and temperatures decreasing towards higher latitudes. Figure B2. Annually averaged difference in solar insolation (left) and temperature (right) caused by extremes in the orbital parameters as listed in Table 5  Obliquity changes produce the largest difference, whereas the precession of the perihelion averages out throughout the year and has no effect on annual averages. Both eccentricity and obliquity show symmetric patterns. However, eccentricity affects the low latitudes the most, meanwhile the difference in the obliquity experiments grows polewards.

Appendix C: Seasonality in the model
The analysis of seasonality discussed in Sec. 3.1 fitted sinusoids to the detrended temperature time series of TransEBM and ERA20CM for each grid box and then extracted the amplitude and phase lag from these fits. The computation of the phase lags relies on the fits, since TransEBM's low temporal resolution would only allow a coarse estimate. However, there are some grid 470 boxes where such a fit oversimplifies the temperature pattern, especially in the ERA20CM reanalysis. This is why the phase lag analysis shown in Fig. 11 excluded some grid boxes. Figure C1 delves deeper into these issues by showing the disagreement between the fitted amplitude and the amplitude calculated as the difference between the minimal and maximal temperature.
The difference between the min-max and fit amplitudes serves as an indicator for how suitable the fit is for deriving amplitudes.
The figure also shows the time series with fits for six exemplary grid boxes.

475
Overall, the difference between the amplitudes is small in TransEBM. The fitted amplitude is larger by up to 1 • over North America and parts of Asia. The min-max amplitude is larger in polar regions, as shown at locations 5 and 6 in Antarctica.
In ERA20CM reanalysis data the differences between the two types of amplitudes are larger, especially in locations with sea ice where they can over 15 • . It is notable that the fitted amplitude is almost always smaller and only ever larger by a little.
This underestimation of the fitted amplitude is visible in locations 5 and 6, where the fit is unable to capture the temperature 480 maxima.
Location 4 near the equator in South America is an example of a grid box that was excluded from the phase lag analysis. It is a location, where the variations have a period smaller than twelve months. Since this was taken as a sign that a sine would not fit the underlying time series, such places where excluded from the analysis. Figure C2 shows the histograms of the phases and periods for both TransEBM and ERA20CM. A period other than around 12 occurs only in few places, mostly around the 485 equator where the time series shows no seasonal cycle.
The phase histograms were used to define cut-offs to distinguish places where sine fitting was suitable from those, where the time series showed little or no sinusoidal behavior. Placing the cut-off is straightforward in TransEBM, since there is a gap at a phase lag of 108 days, after which there are almost no occurrences. Furthermore, Fig. 12a shows that all of these occurrences are close to the equator, that is in places like location 4 without a twelve month seasonal cycle. In ERA20CM, the picture is more complicated as Fig. C2b shows. Here, the cut-off was set by lowering it until all grid boxes close to the equator where the time series is not suited to fitting a sine were excluded. This leads to a cut-off at 150 days. In doing so, large parts of Antarctica were excluded as well. As discussed previously in Sec. 3.1, these are places where the temperature maximum occurs slightly before the maximum in insolation leading to very high computed phase lags, which the analysis omits. While this means that some grid boxes might have been excluded unnecessarily, this procedure allows 495 confidence in all grid boxes that are included.
Author contributions. EZ and KR conceptualized this study, decided on the methodology, created the analyses, wrote, edited and reviewed this manuscript. EZ implemented the software and computations and performed the visualization of the data, which KR supervised.
Competing interests. The authors declare that there are no competing interests.