Evolution of the climate in the next million years: A reduced- complexity model for glacial cycles and impact of anthropogenic CO2 emissions

We propose a reduced-complexity process-based model for the long-term evolution of the global ice volume, atmospheric CO2 concentration and global mean temperature. The model only external forcings are the orbital forcing and anthropogenic CO2 cumulative emissions. The model consists of a system of three coupled non-linear differential equations, 10 representing physical mechanisms relevant for the evolution of the Climate – Ice Sheets – Carbon cycle System in timescales longer than thousands of years. The model is successful in reproducing the glacial-interglacial cycles of the last 800 kyr, in good agreement with the timing and amplitude of paleorecord fluctuations, with the best correlation between modelled and paleo global ice volume of 0.86. Using different model realisations, we produce a probabilistic forecast of the evolution of the Earth system over the next 1 million years under natural and several fossil-fuel CO2 release scenarios. In the natural scenario, 15 the model assigns high probability of occurrence of long interglacials in the periods between present and 120 kyr after present, and between 400 kyr and 500 kyr after present. The next glacial inception is most likely to occur ~ 50 kyr after present with full glacial conditions developing ~ 90 kyr after present. The model shows that even already achieved cumulative CO2 anthropogenic emissions (500 PgC) are capable of affecting climate evolution for up to half million years, indicating that the beginning of the next glaciation is highly unlikely in the next 120 kyr. High cumulative anthropogenic CO2 emissions (3000 20 PgC or higher), which could potentially be achieved in the next two to three centuries if humanity does not curb the usage of fossil-fuels, will most likely provoke Northern Hemisphere landmass ice-free conditions throughout the next half million years, postponing the natural occurrence of the next glacial inception to 600 kyr after present or later.


Introduction
Numerous paleoclimate records show that climate of the last 2.7 Myr is characterized by persistent alternations of glacial and 25 interglacial episodes. The last 800 kyr were dominated by strongly asymmetric glacial cycles with average periodicity of about 100 kyr, a rather slow glacial build-up and much more rapid deglaciations (Hays et al. 1976;Lisiecki and Raymo, 2005).
Antarctic ice core records also show that atmospheric Carbon dioxide (CO2) concentration fluctuated nearly synchronously with the global ice volume, and that CO2 concentration during glacial times was up to 100 ppm lower than during preindustrial time (Petit et al., 1999, Lüthi et al., 2008. 30 The astronomical theory of glacial cycles (Milankovitch, 1941) postulates that growth and shrinkage of ice sheets is primarily controlled by changes in Earth's orbital parameters (eccentricity, obliquity and precession) through their effect on the amount developed in Lord et al. (2016) based on simulations performed with the Earth system model of intermediate complexity cGENIE (Lenton et al., 2006).
The manuscript is organized as follows. In Sec. 2 we introduce the model and datasets. Results are presented in sections 3 and 70 4. Section 3 is devoted to the assessment of model performance. In Sec. 4 we present a probabilistic forecast for the future 1 Myr, under scenarios with and without anthropogenic influence. Finally, in Sec. 5 we summarize the results and present the conclusions.

2
Model and datasets

Modelling approach 75
To develop a probabilistic prediction of future climate we develop a process-based model for the coupled evolution of global ice volume, atmospheric CO2 concentration and global mean surface air temperature. For the past, the model's only external forcing is orbital forcing represented by the maximum summer insolation at 65 o N. For the future, we accounted additionally for the impact of the cumulative anthropogenic CO2 emissions on atmospheric CO2 concentration. The model consists of a system of three coupled, non-linear differential equations, describing Earth system evolution at the timescales from thousands 80 to millions of years. Using different sets of model parameters, we generate a probabilistic forecast of the evolution of the system over the next 1 million years under natural and anthropogenic CO2 release scenarios.
The model contains 12 parameters (see Table 1). Several of them were fixed using prior information, while others were treated as "tunable" model parameters (see Sec. 2.5). For model calibration and testing, we use paleoclimate reconstructions for the 85 last 800 kyr (see below). This period was selected because it is dominated by the long glacial cycles which are expected to continue in the future, at least, for some time (see discussion below) in the absence of significant anthropogenic influence.
Using empirical data and by varying model parameters, we generate a large subset of model realizations of which we select a subset with the highest correlations to paleodata. This subset has been then additionally reduced by requiring consistency with the results of Ganopolski et al. (2016). This reduced subset of model realizations has been finally applied to modelling of the 90 next million years for the natural (no anthropogenic CO2 emission) and several anthropogenic scenarios.

Equation for global ice volume
The first model equation describes the temporal evolution of global ice volume v expressed in nondimensional units (with a value of 1 representative of Last Glacial Maximum (LGM) ice volume levels, according to empirical data). Since it is set that at present v=0, this variable represent anomaly of global ice volume relative to the preindustrial state. As the principal natural 95 forcing for global ice volume is summer NH insolation, it is natural to interpret it as the equation for only NH ice volume. To account for the contribution of the Southern Hemisphere (SH), we use the same approach as in Ganopolski and Calov (2011)  and assume that the ice volume anomaly in the SH closely follows the NH anomaly and makes up a constant, relatively small faction of the global ice volume anomaly. This approach, obviously, is not applicable for a possible future Antarctic and Greenland melting under high CO2 concentrations. This is why we do not consider future sea level rise above the preindustrial 100 level and it is required that v≥0 at any time (see below). The mass-conservation equation states that: "Accumulation" represents here the global ice accumulation (in mass per time units) and "Ablation" represents all mass losses 105 from ice sheets including surface and basal melt, and icebergs calving. We assume that the total accumulation is proportional to the total ice sheet area (Ganopolski et al., 2010). In turn, ice sheet area is closely related to the ice volume. It is often assumed that = , where γ is about 0.8. Here for simplicity we assume γ=1, and thus: where b is the proportionality constant.
On the other hand, the ablation is controlled by several processes of which the energy balance at the surface of the ice sheets is the most important one. Total ablation depends on the size of the ice sheets (and in the NH, especially, on the position of 115 their southern margin), summer insolation (orbital forcing) and CO2 atmospheric concentration (CH4 and N2O also play a role but we assume their radiative forcing to be proportional to the radiative forcing of CO2 as in Willeit et al., 2019). The effect of CO2 on the mass balance of ice sheet is introduced via the logarithm of its concentration because radiative forcing of CO2 is proportional to the logarithm of concentration. As the metric for orbital forcing, as in Ganopolski (2021), we use the maximum summer insolation at 65 o N computed using Laskar et al. (2004). Similar to Ganopolski et al. (2021), we introduced a term 120 depending on dv/dt representing the "domino effect" during glacial termination. To ensure existence of at least one equilibrium solution (solution with dv/dt = 0) for any combination of orbital and CO2 forcing a negative term proportional to the ice volume in power large than one is required (Calov and Ganopolski, 2005;Abe-Ouchi et al. 2013). While Ganopolski (2021) used term v 2 , here we use the term v 3/2 since it gives a better model performance. Thus, the total ablation is represented as: where b01, b2, … b6 are constants, f is the average orbital forcing, and Mv is a memory term that reflects the history of the ice volume during the last yr, as long as dv/dt is negative. At any time t, Mv(t) is defined as: Finally, setting b1=b-b01 and b6 = b06/(1+b5Mv) the mass-balance Eq.

Equation for atmospheric CO2 concentration
It is generally recognized that CO2 (together with other well-mixed greenhouse gases such as CH4 and N2O) represents an 140 important amplifier of glacial cycles and that the anthropogenic CO2 emission will affect the climate during long periods of time (Archer and Brovkin, 2008) and, in particular, the timing and magnitude of future glacial cycles (Archer and Ganopolski, 2005). Although the mechanisms of glacial-interglacial atmospheric CO2 variability are still not fully understood, significant progress in modelling the global Carbon cycle operation during glacial cycles has been achieved in recent times (Brovkin et al., 2012;Menviel et al. 2012;Willeit et al. 2019). Since CO2 and ice volume are highly correlated (during the past 800 kyr 145 the coefficient of correlation between atmospheric CO2 concentration and ice volume is -0.71), it is not surprising that simple conceptual models of glacial cycles can describe the ice volume evolution without an explicit treatment of CO2 (e.g . Paillard 1998;Crucifix 2013;Ganopolski 2021). However, this close relationship between CO2 and global ice volume is not valid for the Anthropocene and, thus, the modelling of the future climate evolution requires an explicit treatment of atmospheric CO2 concentration which is the second equation of our model. The response of the global Carbon cycle to external perturbations 150 involves a broad range of timescales: from very short (annual) to geological. For the purpose of this study we treat natural and anthropogenic anomalies of CO2 separately but consider them additive. Namely, we assume that on the relevant timescales (10 3 -10 5 yrs), the natural component of CO2 concentration is in equilibrium with external conditions and can be expressed through a linear combination of global temperature and global ice volume.

155
The most direct effect of temperature on CO2 is through changes in the CO2 solubility in ocean water (Zeebe and Wolf-Gladrow, 2001). Temperature, directly and through sea ice, also affects ocean circulation (Watson et al., 2015), ventilation rate of the deep water (Kobayashi et al., 2015), changes in relative volume of different water masses (Brovkin et al., 2007) and metabolic rates of living organisms (Eppley, 1972;Laws et al., 2000). The direct effect of ice volume on CO2 changes is less straightforward (please note that the strong effect of ice sheets on atmospheric temperature is already accounted for in the temperature term). An increase in ice volume leads to a decrease of the ocean volume and, as a consequence, increased ocean salinity. As CO2 solubility in the ocean decreases with salinity this would provoke to an increase of atmospheric CO2 concentration, counteracting the temperature effect (Zeebe and Wolf- Gladrow, 2001). Higher global ice volume would also mean a smaller area of the globe is covered by forests and, therefore, a 165 diminished carbon storage in the terrestrial reservoir, driving an increase in atmospheric CO2 (Prentice et al., 2011). Increased glacial supply of iron-rich dust may help suppress the iron limitation in some areas (in particular in the south Atlantic Ocean) leading to increased biological productivity and thus lower atmospheric CO2 levels (Martin, 1990;Watson et al., 2000).
In addition, based on results of Ganopolski and Brovkin (2017) we parameterized CO2 "overshoots" during glacial terminations 170 by an additional term proportional to dv/dt. Note that this additional term is only applied when dv/dt<0. The justification for such parameterization is the results seen in many models (Gottschalk et al., 2019) that a shutdown of the Atlantic Meridional Overturning Circulation (AMOC) causes an atmospheric CO2 rise of about 10-20 ppm. A shutdown of the AMOC is caused by a large meltwater flux which is, in turn, controlled by changes in ice volume (i.e. dv/dt).

175
To describe the anthropogenic component of CO2, we assume that anthropogenic CO2 emission is a relatively short (order of 10 2 yrs) pulse followed by zero emissions. In this case, the temporal dynamics of the anthropogenic component of CO2 anomaly depends only on the cumulative Carbon emission. To describe the long anthropogenic tail of CO2 (AnthCO2 term in the Eq. (7) below) we use the analytical parameterization defined in Lord et al. (2016) based on the results of experiments with the Earth system model of intermediate complexity cGENIE (Lenton et al., 2006). In addition, we assume that natural and 180 anthropogenic CO2 anomalies can be simply summed up and that at the preindustrial time the global Carbon cycle was in equilibrium. This is, obviously, a very strong assumption since even a rather small imbalance in the global Carbon cycle which is impossible to detect at the millennial timescales can result in a very large "drift" of the Earth system from its preindustrial state at the million years timescale. However, due to the absence of any practical alternative, we proceeded with this assumption. 185 After these considerations, the equation that governs the CO2 evolution has the following shape: 190 Where c1… c4 are adjustable model parameters. The first three terms on the right-hand-side of equation (7) represent the already discussed effects of global temperature, ice volume and AMOC weakening or shutdown. The fourth term (c4) is simply a constant and the last term represents external forcing due to anthropogenic processes (and we are, thus, assuming that https://doi.org/10.5194/esd-2021-2 Preprint. anthropogenic emissions occur on top of the natural Carbon cycle as a linear addition). Please, refer to Sec. 4.1 for a detailed explanation on the assumptions and shape of the function AnthCO2(t). 195

Equation for global mean surface temperature
The last equation in our model describes global mean surface air temperature anomaly (∆; relative to preindustrial values) as a linear function of global ice volume and logarithm of CO2 concentration: where d1 and d2 are adjustable model parameters.

Model constraints, bounds and parameter selection strategy
The model of past and future glacial cycles is represented by the system of three equations (6), (7), (8) which contains three variables v, CO2 and ∆T, orbital forcing f(t), anthropogenic CO2 anomaly AnthCO2(t) and twelve parameters (bi i=1:6; cj j=1:4, 205 dk k=1:2). The orbital forcing f(t) depends only on astronomical parameters (eccentricity, precession and obliquity) which are accurately known for the entire period of interest (Laskar et al., 2004). The function AnthCO2(t) is assumed to be null until t = 0 kyr, afterwards the function depends only on time and cumulative CO2 emission as detailed in Sec. 4.1.
Model parameters are selected to yield solutions that are in a good (see discussion below) agreement with the paleoclimatic 210 information over the last 800kyr (see Sec. 2.6 for details on the paleoclimate reconstructions used). To determine some of model parameters and thus reduce the number of "free parameters" we use three conditions based on paleoclimatic data and current estimations of equilibrium climatic sensitivity.
The first constraint is related to preindustrial conditions. We require that if v=0 and ∆T = 0, then CO2 = 278 ppm. This condition 215 fixes the value for c4 at 278.
As the confidence on empirical data for glacial conditions is higher at LGM than at any other glacial episode, the reproduction of LGM conditions is the only constraint that we will require from the model. At LGM, empirical ice volume in nondimensional units was 1, CO2 concentration 194 ppm and global cooling around 5°C (Schneider et al., 2006;Annan et al., 220 2013;Tierney et al., 2020). We then require from the model that: The last imposed criteria is related to the current estimates for the equilibrium climatic sensitivity, i.e. global temperature 225 response to a doubling in atmospheric CO2 concentration from preindustrial levels. We select a climate sensitivity equal to 3.9 C, which coincides with the multimodel mean in the Coupled Model Intercomparison Project 6 (CMIP6; Zelinka et al., 2020): 230 Equations (10) and (11) determine the values for the coefficients d1 and d2 at -3 and 5.56, respectively.

235
Furthermore, to prevent the selection of parameters that will lead to unrealistic values for the selected variables we also apply two limitations.
The first limitation is applied to the simulated ice volume v. First, we assume that for the recent interstadials and any future time, v cannot be negative. Second, it is known that glacial cycles prior to the mid-Brünhes Transition (MBT) about 400 kyr 240 ago differ in some respects from those after MBT. In particular, pre-MBT interglacials were characterised by lower CO2 and higher benthic d 18 O values than the post-MBT episodes. The later implies a large interglacial, probably due to remaining continental ice sheets in the NH. The mechanisms accounting for these differences are still not understood (Tzedakis et al., 2009;Berger et al., 2016) but they are unlikely related to differences in orbital forcing before and after MBT. This is why it is not expected that a simple model forced by orbital variations alone can simulate such behaviour. This is why we prescribe that 245 before the MBT the minimum ice volume must be 0.05 in normalized units: The second limitation, based on empirical data, prevents CO2 to drop to levels significantly lower than the minimum atmospheric CO2 concentration registered in the last 800 kyr by paleorecords (172 ppm) at any given moment in time: 255 In addition, we applied a constraint which is not from the past but from the future. Several modeling studies (Berger and Loutre, 2002;Cochelin et al., 2006;Ganopolski et al., 2016) indicate that the conditions for the new glacial inception will not https://doi.org/10.5194/esd-2021-2 Preprint. Finally, we approach the task of the selection of set of parameters P as a non-linear optimisation problem with equality and inequality constraints. We wish to find P to maximize the optimization target function Cv: 265 where: ( ) � , � denotes the modelled ice volume time-series in the period � , � (i.e. solution obtained for v when using the parameter set P in the system (6-8)); � , � denotes the paleo ice volume time-series in the period � , �; 270 corr(x,y) denotes the linear Person correlation between x and y. The time interval for the optimisation is set to � , � = [−800 , 0 ]. See Appendix A for a discussion of the dependence of model performance on the choice of this time interval. To select parameters that will optimise correlation at the same time as providing magnitudes in accordance to empirical estimations, an equality constraint is enforced: the maximum ice volume must be equal to 1 within a tolerance of 0.15 (in nondimensional units). Finally, the inequality constraint is given by Eq. (14). 275 The time-step for the calculation of the time-series is 1 kyr. We use the interior-point optimisation algorithm under the Matlab environment. Approximately 1000 optimisation routines were performed, each starting from a different combination of the set of nine adjustable parameters (optimisation seeds).

280
Finally, the ensemble of valid solutions, Valid, is defined as the set: The correlation level of 0.7, although arbitrary, guarantees a good fit to the paleo climatic ice volume record. For the selection 285 of solutions, no conditions are imposed on the goodness of fit between modelled and paleo CO2 or temperature. For each P, we calculate v, CO2 and ∆T obtained using P in the system of equations (6), (7), (8) and denote them: ( ) , 2 ( ) , ( ) . This ensemble provides the basis to generate the probabilistic forecast in the remainder of the paper.

2.6
Empirical datasets 290 Paleo reconstructions covering the period [-800 kyr, 0 kyr] are used as part of the learning and/or validation sets for the model.
Reconstructed global ice volume is derived from sea level stack (Spratt and Lisiecki, 2016). Past CO2 atmospheric concentration levels are derived from ice cores records (Lüthi et al., 2008). For global mean surface temperature anomalies (with respect to preindustrial conditions) we use two reconstructions: 1.

Model performance
In this section we analyse the model's performance considering the parameter combinations P belonging to the Valid ensemble.
First, we evaluate the model performance by comparison of model simulations to paleoclimate reconstructions. Second, we analyse the critical insolation -CO2 relationship during glacial inception episodes for the different model realizations derived 300 from Valid and compare them with Ganopolski et al. (2016).
For each P in Valid (353 in total), we calculate ( ) , 2 ( ) , ( ) .  Table 1. By definition, all the solutions derived from the Valid ensemble closely follow the paleo ice volume curve (Fig. 1a), with the ensemble-mean correlation between modelled and paleo ice volume in [-800 kyr, 0 kyr] equal to 0.76. While most of the 310 solutions derived from Valid succeed in capturing the timing and magnitude of the major glaciations, there is a tendency for them to overestimate the ice volume in MIS 18 and MIS 14 (-560 kyr and -520 kyr, respectively). . As model parameters were chosen to maximise ice volume correlation only, it is expected that the performance in terms of CO2 and global temperature is inferior than for global ice volume. The ensemble-mean correlation between modelled and paleo CO2 is 0.5 and some solutions display an amplitude range significantly larger than the observed one, reaching the imposed lower limit of 150 ppm (Fig. 1b). 315

Comparison with paleoclimate reconstructions
For global mean surface temperature (Fig. 1c), the ensemble-mean correlation between modelled and paleo record is 0.56 and the solutions amplitude range within the paleo estimation limits.  While the model demonstrates a good performance when tested against the same data that was used for its training, its ability to produce predictions (situation in which previously unseen data is used as model input) can only be evaluated by the use of independent training and validation samples. A robust estimation of the model's predictive skill, using disjoint training and 325 validation sets, can be found in Appendix A. In general, we conclude that the model has a satisfactory ability also when used in predictive mode and, thus, we confidently venture to utilize it as a tool for the forecast of the next 1 Myr climatic evolution.

Glacial inception: Critical insolation -CO2 relationship
Using a series of CLIMBER-2 model realizations consistent with paleo climatic constraints, Ganopolski et al. (2016) found that the critical threshold for summer insolation at 65° N to trigger glacial inception follows a functional relationship with CO2 345 atmospheric concentration. Given that the radiative forcing of CO2 is proportional to the logarithm of CO2 concentration and that in the CLIMBER-2 the temperature response to CO2 and orbital forcing is approximately linear within the considered range of CO2 concentrations, the critical summer insolation at 65° N can be described as:

=
( 2 280 ⁄ ) + (17) 350 where K = −77 W m -2 and R = 466 W m -2 . Following their approach, we calculate the corresponding coefficient K for each of the ensemble members in the Valid set (Fig. 2). Our results indicate that with the model derived in this study the possible values of the coefficient K range between -1279 and -31 W m -2 , with a median of -393 W m -2 . This highlights that, even though all the solutions derived from the Valid set have a high level of agreement with the paleoclimatic records along last 800 kyr, 355 the relationship between fcritial and log(CO2) is not well constrained by paleoclimate data, i.e model realisations with complitely different values of K can produce similarly good agreement with the paleoclimate reconstructions. However, for the future simulations different K values lead to very different impact of anthropogenic CO2 emissions on the Earth system evolution (see Fig. 3 and the discussion below). This is why we propose carrying forward only those solutions in Valid which are consistent with Ganopolski et al. (2016), i.e. satisfy the condition K ≥ -150 W m -2 . With this choice we tried to account for 360 the fact that this value can be strongly model-dependent. This sub-ensemble, which we call Accepted, consists of 29 elements.
A summary of the characteristics of the Valid and Accepted ensembles can be found in Table 2. The probabilistic distribution of solutions for v, CO2 and ∆T derived from Accepted is displayed in Fig. 4 and the ranges of the different model parameters across Accepted can be found in Table 1.

365
In addition, we define the Best Solution PBest as: PBest values can be find in Table 1. The time-series the whole [-800 kyr, 0 kyr] period. The model skill in reproducing both the temporal variability and amplitudes of the timeseries fluctuations is remarkable. In particular, the performance for ice volume is excellent with a correlation between model and paleodata of 0.86 (Fig. 4a). The model is able to correctly capture the timing of all the glaciations, being the LGM the highest volume event. The ice volume model and paleorecord discrepancy is largest for MIS 18 ([-750 kyr, -710 kyr]). The 375 frequency spectra of ( ) and the empirical estimation are also in good agreement (Fig. 5): there is a clear dominant peak at 100 kyr, with additional power at 41 and 23 kyr. For atmospheric CO2 concentration, the correlation between modelled and recorded series is 0.62 (Fig. 4b). While the model succeeds in reproducing glacial-interglacial variability, it overestimates its magnitude: around 130 ppm instead of 80 ppm. It is, in particular, evident that during interglacials the CO2 restores to its equilibrium value of 278 ppm, which according to the paleorecords is not accurate before the MBT. Finally, the correlation 380 between modelled and paleo global mean surface temperature anomalies is 0.54 (Fig. 4c). The model slightly underestimates the variability of this variable, ranging between 0°C and -5°C, while the observational information suggests a larger range between 2°C and -7°C. The high positive temperature anomalies during some previous interglacials, however, are questionable.

Next 1 Myr
In this section we produce probabilistic forecasts of the evolution of the Earth system in the next 1 Myr, considering natural and anthropogenic emissions scenarios. For the probabilistic forecast we evaluate the future evolution of the solutions derived 395 from the Accepted ensemble.

Description of scenarios
Future scenarios are generated considering different temporal evolution for the excess atmospheric CO2 concentration, which we assume depends only on the cumulative amount of fossil-fuel combustion, AnthCO2(t) (see Eq. (7)). For the natural scenario, AnthCO2(t) is assumed to be zero at all times and, as a consequence, the climate system follows its natural evolution forced 400 only by changes in orbital forcing.   For the fossil-fuel emission scenarios, we consider instantaneous releases of CO2 at t = 0, of different magnitudes. If the duration of the emission pulse release is rather short (centuries) and is followed by zero CO2 emission, the future long-term 415 evolution of CO2 depends primarily on the cumulative emission magnitude while the rate of release is only of secondary relevance (Eby et al., 2009) Fig. 6. 100 kyr after the (low, intermediate or high) emissions pulse release ~ 5% of the initial CO2 anomaly is still present in the atmosphere. For the low emissions scenario, only 500 kyr after the pulse emission the do the remaining anthropogenic CO2 concentration anomalies drop below 1% of the original perturbation, marking the return of the system to its natural evolution. For the medium and high emission scenarios, the return to natural conditions takes even longer. 440

4.2
Forecast under natural scenario Figure 7a displays the natural scenario probabilistic forecast for ice volume. The forecast indicates both periods of high agreement and periods of high discrepancy between solutions, highlighting the co-existence of more certain (deterministic) and lesss certain aspects of the Earth system response to orbital forcing. Most of the solutions agree that the planet will remain in a long interglacial state for the next 50 kyr, during which no large ice-growth is expected. This behaviour is predefined by 450 one of the model constrains, namely, that v should be equal zero at present and by a low orbital forcing during the next 50 kyr.
With high level of confidence, a second long interglacial is expected to occur between 450 and 500 kyr, which is another prolonged period of low orbital forcing related to 400 kyr periodicity of eccentricity. Also with high certainty, short (~10 kyr) interglacials are expected to occur in 110, 230, 250, 270, 370, 620, 820 and 900 kyr from present. Regarding the timing of the next full glacial conditions (defined here as the first moment after present in which ice volume reaches 0.5 in normalized units) 455 the solutions tend to cluster into two different possibilities: predicted to occur in ~90 kyr or in ~150 kyr (see also Fig. 3). For those cases, the ice-growth begins at ~50 kyr or ~110 kyr, respectively. There is also some level of certainty in the timing of occurrence of major glaciations (defined here as periods with ice volume higher than 0.8 in normalized units): in ~100 kyr, ~210 kyr, ~410 kyr, ~450 kyr, ~680 kyr, ~790 kyr, ~870 kyr and 950 kyr from present. The forecast indicates that only two of these major glaciations (the ones in ~210 kyr and ~870 kyr after present) have high probability of reaching LGM magnitudes 460 (~1 in normalized units).
If we consider the forecast produced by the trajectory of the Best Solution (Fig. 8) as the most likely path, under natural circumstances: 1. The next glacial inception is expected to occur ~ 50 kyr after present (with full glacial conditions reached ~90 kyr after present); 2. In the next 1 Myr two long interglacials should take place; 3. Nine major glaciations are expected, 465 all of them with maximum magnitudes around 10-15% lower than LGM level. The natural glacial cycles of the future are strongly dominated by 100 kyr cyclicity as were the past ones.

Forecast under anthropogenic CO2 emissions scenarios
The evolution of the ice volume in the low-emissions anthropogenic scenario (E = 500 PgC) is forecasted to be significantly different from the natural one in the next 200 kyr, indicating a long-lasting impact of fossil-fuel CO2 releases (Fig. 7b). The 470 uncertainty levels for the next 150 kyr and from 400 to 500 kyr are greatly reduced, with most of the solutions agreeing on icefree NH conditions (excepting Greenland). Under this scenario, only after half a million years is the Earth system able to return to an evolution similar to the one expected under natural conditions. Under this setting, the next full glacial conditions are not expected to be able to occur before 180 kyr from present, implying a delay of at least 90 kyr from the natural expected evolution which is consistent with Archer and Ganopolski (2005) and Ganopolski et al. (2016). In particular, the ice volume evolution The probabilistic forecast indicates significantly different conditions from the natural ones in the whole next half million years, with very low probability of considerable ice growth during the totality of this time span (Fig. 7c). In most of the solutions the full glacial conditions occur for the first time either in ~ 550 kyr or between 600 and 700 kyr from present (Fig. 3). In particular, the Best Solution is part of the latter group, with the next full glacial condition taking place only at ~ 670 kyr after present. 485 Finally, with a cumulative pulse emission of 3000 PgC, the Earth system is expected, with high confidence, to experience almost ice-free conditions for the next 600 kyr (Fig. 7d). Under this scenario, the next full glacial conditions are not expected to occur before 670 kyr from present, implying a lag of almost 600 kyr from the natural scenario. For example, in the Best Solution, only two major glaciations are predicted, peaking in 880 kyr and with a magnitude 20% lower than LGM (Figs. 7d 490 and 8a). In the 3000 PgC scenario, half a million years into the future, the CO2 is expected to be 300 ppm (22 ppm higher than preindustrial levels) and the global mean surface temperature still 0.5 o C higher than preindustrial level (Figs. 8b and 8c).

Summary and conclusions
We propose a reduced-complexity process-based model of the evolution of global ice volume, atmospheric CO2 concentration and global mean temperature at multimillennial and orbital timescales. The model only external forcings are the orbital forcing 495 and (for the future) anthropogenic CO2 cumulative emissions. The model consists of a system of three coupled non-linear differential equations representing physical mechanisms relevant for the evolution of the Climate -Ice Sheets -Carbon cycle System in timescales longer than thousands of years. Model parameters are selected to guarantee the best possible agreement with reconstructed ice volume variations during past 800 kyr. Global temperature anomaly at the Last Glacial Maximum and the current best estimate for the equilibrium climate sensitivity have been used as additional constraints on model parameters. 500 Furthermore, the solutions are bounded to account for the difference between pre-and post-mid-Brühnes transition interglacial conditions and a minimum for CO2 concentration is set to avoid unrealistic solutions.
The model is successful in reproducing the natural glacial-interglacial cycles of the last 800 kyr, in agreement with paleorecords, replicating both the timing and amplitude of the major fluctuations, with a correlation between modelled and 505 reconstructed global ice volume of up to 0.86. The model also shows a promising performance when evaluated against data not used for training, indicating a certain skill in predictive mode.
Using different model realizations we generate a probabilistic forecast of the evolution of the system over the next one million years under natural and three different anthropogenic CO2 release scenarios. 510   In the natural scenario, in which anthropogenic CO2 emissions are assumed to be zero at all times, the model assigns high probability of occurrence of: 1. Long interglacials in two periods: [0 kyr, 120 kyr] and [400 kyr, 500 kyr] and 2. The next glacial inception will most likely occur ~50 kyr after present with the development of full glacial conditions ~90 kyr after present. It is also clear, however, that even though there is a high level of agreement in the solutions' trajectories during the 530 past 800 kyr, their paths tend to diverge for the future indicating that the past does not perfectly constraint the future evolution of the climate -ice sheets -Carbon cycle system.
The selected model versions exhibit a large sensitivity to fossil-fuel CO2 releases, with even already achieved emissions (500 PgC) causing a behaviour significantly different from the natural one over extremely long periods of time. The model predicts 535 that a fossil-fuel CO2 release of 500 PgC will make the beginning of the next Ice Age highly unlikely in the next 120 kyr and a delay of half a million years is a possibility that cannot be ruled out. Cumulative fossil-fuel CO2 emissions of the order of 3000 PgC (which could be achieved in the next two to three centuries if humanity does not constraint fossil-fuel usage; Meinshausen et al., 2011) will most likely lead to ice-free conditions over the NH landmasses throughout the next half million years. Under this scenario, the probability for the next glacial inception to occur before 600 kyr is extremely low. Thus our 540 study demonstrates that, in spite of large unecrtainties in model parameters, there is a qualitative difference in the long-term Earth system trajectories (Steffen et al., 2018) between cumulative emissions of 500 and 3000 GtC.
The present study reveals one serious challenge in modeling of the deep future: the results show very strong sensitivity to the relationship between critical insolation threshold for glacial inception and CO2 levels. In turn, this relationship is poorly 545 constrained by the paleoclimatic data because during previous interglacials CO2 was close to or lower than the preindustrial level. So far, this relationship has been analysed only with the CLIMBER-2 model (Calov and Ganopolski, 2005;Ganopolski et al., 2016). This is why in this work, we assumed the uncertainty in this relationship of 100%. Reducing this uncertainty by performing experiments similar to those described in Ganopolski et al. (2016) but with more advanced Earth system models can help to reduce uncertainties in future projections. 550

Appendix A: Estimation of predictive skill
When a model is developed with a predictive goal (like in this study), it is important to have an assessment of its predictive 555 skill. The best possible estimations of this skill are obtained by evaluating how well the model performs when presented with previously unseen data.
To follow the cross-validation assessment criterion (Stone, 1974), we divide the period [-800 kyr, 0 kyr] in halves: P1 = [-800 ky, -400 kyr] and P2 = [-400 kyr, 0 kyr]. First, P1 is labelled training set and P2 validation set. Second, the model is trained 560 using the information in P1: we run the Optimisation scheme to maximise the correlation between modelled and paleo ice volume during P1; we then accept all the models (parameter combinations) yielding correlations higher than 0.7. Third, all the accepted models are validated using the information in P2. We select three different performance metrics: the correlation between modelled and paleo ice volume, between modelled and paleo CO2 or between modelled and global mean surface temperature anomaly, calculated during P2. Finally, we repeat the previous three steps but labelling P1 as validation and P2 as 565 training sets, respectively. For each metric, the average that the models achieve on the two possible validation sets considered is the estimate of the model's predictive skill.
In Table A1 we summarize the results of this cross-validation assessment, showing the results for the three performance metrics. For reference, Table A1 also displays the values of the metrics when calculated over the training period (instead over 570 the validation period). Results are the average of two possible cases: P1 as training and P2 as validation sets, or P2 as training and P1 as validation sets. As expected, when the same period is used both for model training and evaluation the results are in better agreement with paleoclimate data: the ice volume metric is 0.78, the CO2-metric 0.52 and the temperature-metric 0.58.
These values drop to 0.49, 0.36 and 0.34 respectively, when the independent validation set is used to calculate the metrics.
Even though the drop in performance is not negligible, the model skill is still substantial in particular for predicting ice volume. 575 Table A1: Estimation of Model Predictive Skill. Only solutions leading to a correlation between modelled and paleo estimation of ice volume higher than 0.7, calculated over the training period, are considered. Please see text for more details on the training and validation periods used.

Author contribution 585
AG conceived the original idea. ST performed all the calculations and produced figures. Both authors contributed to the analysis of results and writing of the manuscript.