Coupled atmospHere Interior modeL Intercomparison (CHILI). I. Evolutionary Modelling –
Primordial Magma Oceans of Earth and Venus

This part of CHILI uses Earth and Venus as case studies for comparing simulated interior-atmosphere evolution pathways. We consider ‘nominal’ scenarios for both planets, initialised with $4.7\times 10^{20}$ kg of hydrogen (equivalent to three Earth oceans³³3Here, ‘oceans’ of hydrogen content is a unit describing an amount of hydrogen atoms, equivalent to the mass of hydrogen atoms within Earth’s surface water ocean reservoir. One ocean corresponds to $1.4\text{\times}{10}^{21}\text{\,}\mathrm{kg}$ of \chH2O, and $1.6\text{\times}{10}^{20}\text{\,}\mathrm{kg}$ of hydrogen atoms (Clark, 1982).) and $2.73\times 10^{20}$ kg of carbon in the mantle-atmosphere system. If outgassed entirely as \chCO2, this carbon inventory corresponds to 192.4 bar of surface pressure. Nominal Venus is identical to Nominal Earth except its mass, radius and orbital distance are adjusted to match Venus’ present day values (Lodders & Fegley, 1998).

Magma ocean cooling timescales are sensitive to planetary volatile budgets, particularly that of hydrogen, because thick atmospheres induce strong blanketing and greenhouse effects which delay mantle solidification (Nicholls et al., 2024; Nikolaou et al., 2019; Hamano et al., 2015). Yet, multiple physical processes and circumstantial factors collectively determine the volatiles obtained by a planet during its formation (Ikoma & Genda, 2006; Krijt et al., 2022). For example, a planet could form from volatile-richer or volatile-poorer material, depending on its location exterior or interior to the water ice line or soot line (Drążkowska et al., 2023; Bergin et al., 2026), while being subject to desiccation and mixing of multiple protoplanetary disk sources (Lichtenberg et al., 2019; Ikoma & Genda, 2006; Sossi et al., 2025). Giant impacts, later-stage delivery, and volatile redistribution would then modulate the amount of volatile elements distributed across early mantles and atmospheres (Lupu et al., 2014; Roche et al., 2025; Nakajima et al., 2021; Kegerreis et al., 2020; Lock & Stewart, 2024). As a result, there are large empirical uncertainties on Earth’s initial (Wang et al., 2018) and present volatile budget: its combined mantle and atmosphere water content is estimated between $\sim 2$ and $\sim 10$ oceans (Peslier et al., 2017). There is substantial uncertainty in the additional amounts of species of light elements incorporated in Earth’s metallic iron core (Hirose et al., 2013; Shahar et al., 2026). To respect the uncertainty and potential range of outcomes arising from these processes, we additionally consider a range of evolution scenarios, across three hydrogen and three carbon inventories: [1.60, 7.80, 16.00] $\times~10^{20}$ kg H and [1.36, 2.73, 5.44] $\times~10^{20}$ kg C, respectively. These hydrogen budgets correspond to 1, 5, and 10 oceans of \chH2O (Peslier et al., 2017; Javoy, 2005). These elemental budgets represent the light elements in the mantle plus atmospheres that are free to cycle between these two reservoirs. No chemical interaction between the iron core and planetary mantle is taken into account here. Our exploration of various initial hydrogen budgets reflects our uncertainty on the amount of \chH2O and \chH2 to be exogenously delivered (Ormel et al., 2021; Krijt et al., 2022) and endogenously produced (Ikoma & Genda, 2006; Genda, 2016) within young planets. We anticipate a range of behaviours to arise from these nine cases, since the CHILI models incorporate various physical parametrisations of atmospheric opacity, chemistry, and partitioning (Lichtenberg et al., 2026).

All simulations treat the metallic core as chemically inert, with a fixed 55% interior radius fraction informed by Earth’s interior structure (Lodders & Fegley, 1998). We initialise simulations at a stellar age of 50 Myr (Rubie et al., 2025; Dalrymple, 2001). The models adopt different initial temperatures and entropies, while ensuring that these correspond to a fully molten state. This initial condition ensures that the planet’s core has fully differentiated from the mantle, the Sun has entered its main sequence phase, and the solar nebula has dispersed (Halliday & Canup, 2023; Schoenberg et al., 2002; Baraffe et al., 2015). Additional accretion of material onto these young planets is not simulated by any of the CHILI models, but would be minimal after 50 Myr (Dauphas, 2017; Halliday & Canup, 2023).

To curtail model complexity, all of these simulations neglect clouds or aerosols and their radiative effects. To enable intuitive model intercomparison, they also disable Rayleigh scattering in their radiative transfer calculations, and instead enforce a constant 10% Bond albedo scale factor applied to incoming stellar radiation (Pluriel et al., 2019; Cmiel et al., 2025). Various treatments of atmospheric escape are applied, so a 30% energy-limited photoevaporative escape efficiency is adopted where applicable. Simulations terminate according to model-specific criteria on mantle melt fraction.

The CHILI models’ configuration files and the resultant simulation data, in addition to intercomparison analysis scripts and plots, are freely available on CHILI’s GitHub repository⁴⁴4https://github.com/projectcuisines/chili/. This GitHub repository is archived on Zenodo⁵⁵5https://zenodo.org/records/20680020: Nicholls (2026a).

The eight evolutionary models compared in this paper are outlined by Table 1. This suite differs from the CHILI protocol by three key aspects: (i) we consider only explicit time-tracking codes, since static models are left for CHILI Part III; (ii) with the addition of NEONGOOEY as a successor to GOOEY; (iii) with the inclusion of the MOAI, LINCS, and PlanAtMO as additional codes. The following paragraphs briefly outline the important aspects of the models participating in this paper’s intercomparison and highlights their key differences.

Figure 1 compares nominal Earth and Venus evolution tracks of volumetric mantle melt fraction, calculated by each model (panels). The time $x$-axis represents the simulated time from an initially molten state, such as following a particular giant impact event (Abe, 1997). All of these Nominal Earth cases qualitatively agree by cooling and terminating (solid scatter points) before two empirical estimates of maximum magma ocean solidification time from surviving Hadean rocks and zircons (vertical dashed lines; Sole et al. (2025)), with NEONGOOEY taking the longest time to solidify (3.8 Myr).

Venus cases (dashed lines) typically have longer solidification timescales than Earth; its enhanced irradiation flux offsets thermal emission to inhibit planetary cooling. However, the Earth-Venus evolutionary divergences vary between the models (solid vs dashed lines). Simulations of Nominal Venus, except by MOAI and GOOEY, enter into states of global radiative equilibrium in which stellar radiation absorption balances thermal emission, and the melt fraction becomes approximately constant in time (van Dijk et al., 2026; Nicholls et al., 2024). We may classify these as ‘Type II’ planet behaviours in contrast to the ‘Type I’ scenarios which continually cool to the point of solidification (Hamano et al., 2013). During these periods of radiative equilibrium, mantle melt fraction is approximately constant and the simulated Venusian primordial magma oceans are extended for up to 49 Myr (dashed green line). Type II scenarios eventually solidify after long-term atmospheric escape processes sufficiently strip the atmospheric greenhouse. LINCS and PROTEUS (green and orange lines) find markedly different evolutionary pathways between their Nominal Venus and Nominal Earth scenarios: predicting Type I behaviour with early solidification for Earth, and Type II behaviour with late solidification for Venus. Magma oceans facilitate efficient interior-atmosphere exchange of volatiles, so Type II scenarios with extended magma ocean periods could offer greater exposure of chemical (e.g. H content) and isotopic signatures (e.g. D/H) to atmospheric escape processes (Grinspoon, 1993). The modelling approaches which give rise to these different outcomes are identified in later sections.

Our comparisons between Nominal Earth and Venus show a range of divergent thermal evolution pathways for these two planets, depending on the model applied to simulate them. Since post-formation volatile budgets are expected to vary between planets (Krijt et al., 2022; Drążkowska et al., 2023), and chemistry is suggested to exert a leading-order control over planetary cooling timescales (Nicholls et al., 2024), we investigated how these models’ respond to our grid of initial H and C inventories. Figure 2 plots the simulated time taken for our grid of Earth scenarios to cool to three distinct melt-fraction thresholds: 95% (a, mostly molten), 40% (b, rheological transition), and 5% (c, mostly solid). Circular scatter points show different hydrogen and carbon inventories, across our grid of parameters, compared to the Nominal Earth scenario (crosses).

For a given model (scatter colour), thermal evolution timescales correlate with the hydrogen inventory ($y$-axis): Hhigh cases take longer to reach all melt fraction thresholds (panels) compared to Hmid and Hlow scenarios. This known correlation between hydrogen content and cooling rate has been well established in the literature since additional volatiles (i) introduce a stronger atmospheric blanketing effect (Nicholls et al., 2024; Hamano et al., 2015), and (ii) are sustained against escape processes for a longer time (Hamano et al., 2013; Schaefer et al., 2016). Variations in carbon inventory (marker opacity) have a secondary effect to the hydrogen inventory, which is also an established trend (Nicholls et al., 2024; Miyazaki & Korenaga, 2019), and explained by \chH2O and high-pressure \chH2 having a stronger greenhouse potential than \chCO2 and \chCO (Lichtenberg et al., 2021; Wordsworth & Pierrehumbert, 2013; Ramirez et al., 2014).

The simulated time required for Earth to reach a given melt fraction (Figure 2 panels) is more sensitive to the model applied (scatter colour) than to the volatile inventory considered (scatter opacity). An important corollary being that any inferred thermal histories (e.g. from the estimated age of ancient mineral samples) from initial volatile budgets is currently challenging, because it is more strongly dependent on model selection than on the inferred parameters.

Non-linear cooling sequences are common across the models (Figure 2). Most models initially cool rapidly, with PROTEUS being the slowest to reach 95% melt fraction (panel a), although this still occurs within 0.1 Myr of simulated time. PROTEUS’ departure from the other models, in this manner, is explained by its mantle dynamics treatment (Section III.6) and mantle melting curves (Appendix A). The time taken to reach 40% melt fraction is also rapid on geologic timescales: less than 0.6 Myr across all cases (panel b). However, the thermal evolution stalls at these partially-molten states near the rheological transition between 30–40% melt fraction, across all models (Costa et al., 2009; Abe & Matsui, 1985). This commonality arises because partially molten ‘mushy’ magma ocean layers have increased viscosities with lower convective heat fluxes (Schaefer et al., 2016).

In the following sections, we use diagnostic variables derived from Nominal Earth and Venus simulations to explain these different cooling timescales, which arise from the key physics that shapes these early stages of planetary thermo-compositional evolution.

Cooling timescales and the emergence of Type I/II behaviours are sensitive to atmospheric mass and composition, because energy released from the planet’s interior is largely blanketed by infrared-opaque atmospheric species (Nicholls et al., 2024; Lichtenberg et al., 2021; Zahnle et al., 2010). The diversity in inter-model cooling timescales shown in Figures 1 and 2 potentially reflects their range of calculated atmospheric compositions. To test this, Figure 3 plots outgassed atmospheric compositions (bars) and surface temperatures (stars) for simulations of Nominal Earth during its early molten state (panel a) and nearly solidified state (panel b).

Most models predict early \chCO-\chCO2 dominated atmospheres (yellow and orange bars, panel a) because early molten silicate readily dissolves hydrogen bearing \chOH- ions, while outgassing C-rich atmospheres (Sossi et al., 2023; Elkins-Tanton, 2008; Bower et al., 2022). The total surface pressures of C-rich models (NEONGOOEY, PROTEUS, PACMAN, MOAI, PlanAtMO) show good quantitative agreement, with variation in total surface pressure from $\sim 120$ to $150$ bar. For Nominal Earth at volumetric melt fraction $\Phi=95\%$ (Figure 3a), PROTEUS has the lowest surface temperature $T_{\mathrm{surf}}=2266$ K and produces an almost-entirely \chCO2 atmosphere (88%), while MOAI and PACMAN favour \chCO-dominated compositions. These three different \chCO2/\chCO ratios are partially explained by the models’ different surface temperatures (grey stars) at 95% mantle melt fraction; \ch2 CO + O_2 $\Leftrightarrow$ 2 CO2 is an exothermic reaction which favours \chCO at higher temperatures (Hirschmann, 2012). However, differences in modelled carbon speciation are also driven by the upper-mantle redox state which sets the \chO2 partial pressure (Figure 6). PROTEUS pins $f$\chO2 relative to the iron-wüstite buffer, leading to $26\times$ the $p$\chO2 calculated by PACMAN (at 2750 K, where the mantle is entirely molten), thereby favouring \chCO2 over \chCO in Figure 3a. MOAI predicts a substantially more reduced initial magma ocean state (Figure 6), leading to its predicted $p$\chO2 being only $\sim 1/420$ that of PACMAN’s, explaining MOAI’s prediction of a \chCO-dominated outgassed atmospheric composition (Figure 3a) despite its lower surface temperature than PACMAN. Together, these differences in carbon thermochemistry highlight the intimate relationship between mantle redox state and atmospheric composition.

GOOEY and LINCS outgas thin steam atmospheres (12 and 21 bar, respectively) because they neglect carbon chemistry. NEONGOOEY does not model \chCO but does predict 87% \chCO2-dominated compositions comparable with PROTEUS and PlanAtMO. Since the Nominal Earth atmospheres in GOOEY and LINCS are thinner, they are expected to less effectively ‘blanket’ the planet’s interior, resulting in correspondingly larger outgoing radiation fluxes and shorter solidification times (Figure 1). However, the effect of including the opacity of carbon-bearing species on the planet’s thermal evolution remains a second order effect, compared to the other modelling differences which dominate at these early evolutionary stages (Section III.6). The outgoing radiation flux from GOOEY is comparable with PACMAN and PROTEUS, despite these latter two codes also simulating carbon chemistry and the additional opacities of \chCO2 plus \chCO (Figure 8b). This further demonstrates that inter-model differences yield larger changes to cooling timescale than variations in initial whole-planet carbon budgets (Figure 2).

Atmospheres overlying nearly-solidified mantles (Figure 3b) are generally found to be dominated by steam, but show substantially increased diversity in composition and total pressure, compared to initial fully-molten states. Bulk carbon inventories are largely outgassed at 95% melt fraction, so magma ocean solidification towards 5% melt fraction injects hydrogen atoms into initially C-rich atmospheres (Bower et al., 2022; Nicholls et al., 2024; Nikolaou et al., 2019). At the same time, cooler surface temperatures favour thermochemical conversion of \chCO into \chCO2, resulting in \chH2O+\chCO2 dominated mixtures at later stages (Figures 3b and 5).

Using our simulations of Nominal Earth atmospheric compositions, shown in Figure 3, we are able to identify three classes of behaviours demonstrated by these models’ evolving atmospheric compositions.

• •

  NEONGOOEY, PROTEUS, and PlanAtMO predict early \chCO2-rich atmospheres, with total pressures of approximately $\sim 150$ bar, that become enriched \chH2O which dominates the composition at later magma ocean stages.

• •

  PACMAN and MOAI predict early $\sim 150$ bar CO-rich compositions, that transition to \chCO2-dominated with only minor \chH2O degassing.

• •

  GOOEY and LINCS do not model carbon chemistry, so predict very thin \chH2O atmospheres at early stages, which thicken during magma ocean solidification due to ‘catastrophic’ degassing of otherwise-soluble \chH2O.

These three classes will help us interpret the outputs from these models’ climate calculations, in Section III.5.

At $\Phi=5\%$ for Nominal Earth, PACMAN predicts the lowest surface pressure (176 bar, 71% \chCO2) while PROTEUS predicts the largest (578 bar, 82% \chH2O). This range primarily results from these models’ differential treatment of hydrogen partitioning within the planet’s interior (Carone et al., 2025). To demonstrate this, Figure 4 plots the absolute mass inventories of hydrogen (green) and carbon (orange) distributed by each model between the atmosphere (panel a), molten magma ocean (panel b), and solid mantle regions (panel c). The PACMAN Nominal Earth simulation stores $\approx 50\%$ of its hydrogen inventory in the solidified mantle, through its parametrised dynamic interstitial melt trapping (Laurent et al., 2020; Schiano, 2003), which is theorised to be an efficient mechanism for deep sequestration of hydrogen (Sim et al., 2024; Kent, 2008; Hier-Majumder & Hirschmann, 2017). Similarly, MOAI’s Nominal Earth simulation incorporates 79% of its remaining hydrogen into the deep mantle via the interstitial melt pockets, with minimal (6%) atmospheric storage after substantial mantle solidification and atmospheric escape (Figure 7). These differences are an important source of model disagreement because PACMAN’s and MOAI’s predicted deep hydrogen storage leaves a hydrogen-poor atmosphere (Figure 3b), so H atoms are more readily retained against escape processes. In contrast, PROTEUS entirely neglects volatile storage within the solid phase of the mantle, so PROTEUS predicts massive hydrogen-enhanced atmospheres containing \chH2O partial pressures exceeding 450 bar (Figure 3b), alongside trace amounts of \chH2, at these later stages of magma ocean evolution.

Rather than considering trapped pockets of melt within the solidified mantle regions, all other CHILI models simulate the equilibrium partitioning of \chH2O between the molten and solidified regions of the mantle – where \chH2O molecules are considered to be directly incorporated within the crystalline lattice – parametrised using equilibrium partition coefficients (Lebrun et al., 2013; Elkins-Tanton, 2008). GOOEY, NEONGOOEY, and LINCS adopt the same coefficient value $k_{\mathrm{H_{2}O}}=10^{-2}\mathrm{\,kg/kg}$, which leads to small amounts of \chH2O storage in the solid mantle phases (Figure 4c). The partition coefficient simulated by PlanAtMO depends on the pressure at which solids form, which determines their mineral phase (Schaefer et al., 2024; Walbecq et al., 2025), generally resulting in $k_{\mathrm{H_{2}O}}$ values less than $10^{-3}$. These models show good agreement on the amount of \chH2O sequestered by equilibrium partitioning (Figure 4b,c).

MOAI’s predicted $3\text{\times}{10}^{20}\text{\,}\mathrm{kg}$ of solid-mantle hydrogen storage is consistent with a lower estimate of $3.4\text{\times}{10}^{20}\text{\,}\mathrm{kg}$ of hydrogen atoms (stored as water) within Earth’s present mantle (Peslier et al., 2017). The other models retain substantial hydrogen inventories, comparable to Earth’s present water content, within their remnant magma ocean regions (Figure 4b). Whether these hydrogen reservoirs would be finally outgassed during magma ocean crystallisation, or trapped within the mantle, is not resolved here (Mojzsis et al., 2001; Foley, 2015; O’Neill et al., 2007).

Earth and Venus are similarly sized, and their Nominal scenarios are initialised with the same carbon and hydrogen budgets (Section II). Thus, modelled differences in their thermal and compositional histories are reflective of their differing bolometric and high-energy instellation fluxes. Figure 5 plots the outgassed atmospheric composition and surface temperature of Nominal Venus at $\Phi=5\%$. This scenario is directly analogous to the Nominal Earth case in Figure 3b and broadly exhibits quantitatively similar compositions, for a given model. Comparison of Figures 3b and 5, which represent late magma ocean stages ($\Phi=5\%$), shows that total surface pressures on Nominal Venus are smaller than for Nominal Earth, despite the simulations being initialised with the same carbon and hydrogen budgets. The lower surface pressure on Venus arises because the planet is exposed to larger X-ray and ultraviolet radiation fluxes that drive atmospheric escape (Zahnle & Kasting, 1986; Hunten et al., 1987; Owen, 2019). Since Venus is also exposed to a larger bolometric stellar flux, a thinner atmosphere is necessary for it to solidify (Hamano et al., 2015; Krissansen-Totton et al., 2021b; Hamano et al., 2025).

An exception to these Earth-Venus pressure differences is MOAI, which only considers the escape of \chH2, of which only trace amounts are formed, so MOAI predicts similar surface pressures between Venus and Earth cases (218 and 215 bar, respectively). PACMAN also does not follow this trend of atmospheric pressure difference between Earth and Venus. The PACMAN model finds a surface pressure of 182 bar for Venus compared to 176 bar for Earth at the same melt fraction. This can be attributed to Venus’ smaller mantle volume, and the treatment of melt trapping of volatiles and atmospheric escape in the PACMAN model. Competition between atmospheric radiative blanketing and heating by stellar radiation is consistent with later solidification under the Type II scenario where escape regulates the cooling timescale. Type II Venus is exhibited by all models except MOAI and PACMAN, for which both its Nominal Venus and Nominal Earth scenarios directly cool to solidification (Figure 1).

Hydrodynamic outflow of escaping gas, driven by the absorption of stellar radiation, is expected to preferentially remove lighter atoms under most irradiation regimes (Hunten et al., 1987; Zahnle & Kasting, 1986; Owen, 2019). This would increase the O/H and C/H mixing ratio of the atmosphere (Cherubim et al., 2025; Krissansen-Totton et al., 2021a), and potentially drive mantle oxidation by hydrogen loss (Wordsworth & Pierrehumbert, 2014; Katyal et al., 2020). GOOEY and PACMAN incorporate these effects but show close chemical correspondence between Nominal Venus and Nominal Earth at 5% melt fraction. For example, at $\Phi=5\%$, PACMAN calculates atmospheric C/H mass ratios of 31.25 and 25.46 for Nominal Earth and Venus respectively – indicative of systematic carbon enhancement on Earth, opposite to our expectations from fractionation by atmospheric escape. So, while fractionation by escape is known to imprint measurable isotopic signatures on these planets’ atmospheres, it plays only a minor role in shaping their thermal histories (Zahnle & Kasting, 1986; Grinspoon, 1993). These C/H differences arise through melt trapping within the solid phase of the mantle, as modelled by PACMAN, where faster solidification (in the Earth case) more efficiently traps pockets of melt in the solid phase as the magma ocean solidifies upwards from the core-mantle boundary. This trapping leads to higher atmospheric C/H ratios, given a fixed initial inventory C and H atoms (Carone et al., 2025).

Mantle oxidation state is expected to evolve due to multiple physical effects which can compound and negate each other (Schaefer & Fegley, 2017; Lichtenberg, 2021; Hirschmann, 2022). We may proxy the mantle oxidation state using thermochemical redox reactions; e.g. the iron-wustite buffer \ch2 Fe + O2 $\Leftrightarrow$ 2 FeO. This reaction is important and often adopted as a proxy because (i) iron atoms can take on different oxidation states, (ii) iron is abundant in planets due to its pile-up during stellar nucleosynthesis (Lodders & Fegley, 1998; Kippenhahn et al., 2012), and (iii) oxygen is the most abundant element in rocky planets (Wang et al., 2018; Sossi et al., 2025). For these reasons, the fugacity of oxygen in the upper-mantle ($f$\chO2) can be used as a proxy for the mantle redox state, by its quantification relative to the $f$\chO2 value of the iron-wüstite buffer (Frost, 1991). This reaction is temperature dependent, so iron redox thermochemistry and differential melt-solid partitioning of iron (e.g. between ferric and ferrous states) may sequentially oxidise cooling magma oceans (Kress & Carmichael, 1991; Schaefer et al., 2024). However, this temperature effect is modulated by atmospheric escape of oxygen atoms, mantle dynamics, and core formation, which dynamically change the availability of iron and oxygen atoms in the mantle (Wade & Wood, 2005; Wordsworth et al., 2018; Kasting et al., 1993). Critically, $f$\chO2 also controls chemical speciation in these models’ outgassed atmospheric compositions because it is also treated as the partial pressure of \chO2 molecules; e.g. via \ch2 H2 + O2 $\Leftrightarrow$ 2 H2O and \ch2 CO + O2 $\Leftrightarrow$ 2 CO2.

Figure 6 plots the mantle $f$\chO2 calculated by four simulations of Nominal Venus. In absolute terms, the simulated $f$\chO2 decreases over time as temperature decreases (Figure 6a). However, it is more appropriate to compare $f$\chO2 in relative terms (Frost, 1991). Figure 6b shows four models reproducing the CHILI protocol’s requirement that endpoint $f$\chO2 values tend towards IW+4. However, the models show different temperature-dependent behaviours. PROTEUS (orange line) prescribes $f$\chO2 relative to IW, so its $f$\chO2 values necessarily sit near IW+4 across all temperatures (and corresponding planet ages) and do not respond to O/H fractionation from atmospheric escape, nor to ferric-ferrous iron partitioning during magma ocean crystallisation. PACMAN tracks oxygen mass-conservation, while also accounting for hydrogen-oxygen fractionation in the escaping hydrodynamic outflow. Despite these additional physics and Venus’ enhanced exposure to X-ray and ultraviolet radiation compared to Earth, PACMAN $f$\chO2 values remain near IW+4 across all temperatures: the blue line in Figure 6b shows only slight oxidation from IW+2.8 to IW+3.2 during simulated cooling of the surface temperature from 3250 to 2500 K. Hydrodynamic fractionation of O/H in the PACMAN case is minimal and therefore does not drive mantle oxidation; consistent with PACMAN’s predictions of minimal atmospheric loss (Figure 7), and with its predicted preferential sequestration of hydrogen into the solid-phase of the mantle (Figure 4c), and production of carbon-rich atmospheres (Figure 5).

The oxygen fugacity calculated by GOOEY (grey line in Figure 6) only shifts from IW+0 to IW+4 during the final stages of the simulated cooling of Nominal Venus, at melt fractions less than 5%. This late-stage oxidation is also not driven by hydrodynamic fractionation of O/H – since GOOEY predicts minimal hydrogen loss (Figure 7) – and is caused entirely by the concentration of ferric iron in the remnant magma ocean layers.

In comparison to GOOEY and PACMAN, MOAI shows a more continuous trend of f\chO2 as a function of temperature (and time) – plotted in purple in Figure 6. MOAI’s calculated $f$\chO2 also span the widest range: from initially IW+0 to IW+4 at 5% mantle melt fraction. This behaviour arises, firstly, from MOAI’s prediction of continuous and substantial hydrogen loss (as \chH2), since its simulations retain their entire oxygen budget, which drives monotonic oxidation of the planet’s mantle (Kasting et al., 1993). MOAI also simulates differential ferric-ferrous iron partitioning into the solidified mantle regions, causing oxidation of the magma ocean under these bottom-up fractional crystallisation scenarios (Maurice et al., 2023).

Overall, these differences in volatile chemistry and partitioning highlight that: (a) interiors and atmospheres of Earth- and Venus-like scenarios are predicted to oxidise over time, under bottom-up fractional crystallisation scenarios; (b) processes which sequester hydrogen into deep planetary interiors strongly modulate atmospheric compositions, controlling whether carbon- or water-dominated atmospheres are produced at magma ocean solidification; (c) mantle redox evolution pathways regulate outgassed atmospheric \chCO2/\chCO ratios.

Atmospheric escape incorporates multiple physical, chemical, and hydrodynamic processes, which resolve towards different escape ‘regimes’ depending on a planet’s size and irradiation exposure (Owen, 2019). No specific protocol is established for atmospheric escape; the CHILI models take differing approaches to parametrising atmospheric escape, although they generally adopt the energy-limited approximation in which upper-atmosphere heating by absorption of X-ray and ultraviolet radiation launches a hydrodynamic outflow of gas (Hunten et al., 1987; Zahnle & Kasting, 1986; Owen, 2019). PROTEUS treats the outflow as intrinsically non-fractionating, while MOAI and LINCS assume only \chH2 can escape, PACMAN accounts for C-H-O fractionation by thermospheric diffusion, and GOOEY/NEONGOOEY both treat O-H fractionation. PlanAtMO does not simulate atmospheric escape.

Figure 7 plots the hydrogen and carbon inventories of the Nominal Venus scenario, over time, which both decrease due to the cumulative effects of atmospheric escape processes. The circular markers show the amount of initial H and C retained by the planet by the time its mantle has reached a simulated melt fraction of 5% – corresponding to the outgassed atmospheric compositions shown in Figure 5. LINCS predicts that Venus retained only 22.9% of its post-accretion hydrogen budget at the point of magma ocean solidification (green line), in contrast to GOOEY, which predicts 99.8% hydrogen retention, despite these models both considering pure-\chH2O atmospheres and equilibrium partitioning between solid-mantle, magma ocean, and atmosphere. Their different predictions arise from the prolonged magma ocean in the LINCS case, where Venus exhibits type II behaviour that exposes the planet’s hydrogen content to atmospheric escape for 40 Myr, while GOOEY solidifies in a much shorter timescale of 0.45 Myr without entering a state of global radiative equilibrium. PROTEUS predicts a similar magma ocean cooling timescale to LINCS, while finding that a larger fraction (63.2%) of the initial hydrogen budget is retained. This difference between LINCS and PROTEUS arises because PROTEUS calculates evolving bulk mass-loss escape rates, which are then partitioned into carbon and hydrogen mass fluxes, so PROTEUS’ substantial outgassed carbon reservoirs offset the amount of hydrogen lost. The massive loss of carbon in the PROTEUS case corresponds with its prediction of a steam-dominated (i.e. oxygen and hydrogen) atmospheric composition, shown in Figure 5.

MOAI finds that 20% of the initial hydrogen is removed by escape during the course of simulated Nominal Venus magma ocean evolution, which corresponds with gradual mantle oxidation in Figure 6. PACMAN and GOOEY predict negligible amounts of hydrogen escape in the Nominal Venus case, and correspondingly predict only minor mantle oxidation until the mantle melt fraction reaches 5%, at which point other processes drive substantial changes in $f$\chO2 (Section III.3). Together, these three models’ behaviours suggest that substantial escape of hydrogen does correlate with oxidation of rocky planet interiors, although it acts simultaneously with interior mantle processes that can drive a final, rapid oxidation (Figure 6).

Substantially different predictions of hydrogen and carbon losses between these models, shown by Figure 7, is an important point of model disagreement. Venus’ and Earths’ surficial volatile inventories after magma ocean solidification likely played a key role in controlling their early climate states (Constantinou et al., 2026; Sossi et al., 2020). However, these differences are not unexpected, because accurate treatment of atmospheric escape processes requires vertically resolved hydrodynamic photochemical modelling of these planets’ upper-atmospheres (Schulik & Booth, 2023; Brain et al., 2016; Owen, 2019).

The rate of energy radiated to space by a planet determines its cooling rate. This net energy flux is the balance of outgoing thermal radiation, incoming stellar radiation, and scattering processes (Pierrehumbert, 2010). Thus, we must consider how the CHILI models determine these flux terms in order to understand their predicted planetary cooling timescales (Lichtenberg et al., 2021; Lebrun et al., 2013). Outgoing longwave radiation (OLR) quantifies bolometric radiation emitted by the planet to space, which is determined by thermal emission and attenuation by the atmosphere. The OLR is effectively controlled by the radiating temperature of the planet’s photosphere, where the optical depth of longwave radiation $\tau_{\mathrm{LW}}=1$, the location of which is determined by the atmospheric composition and temperature profiles (Pierrehumbert, 2010; Stamnes et al., 2017). Aerosol scattering effects are not modelled, per the CHILI protocol, so a positive difference between the OLR and absorbed stellar radiation (ASR) represents net energy loss.

Figure 8 plots outgoing longwave radiation fluxes from Nominal Earth simulations, against the modelled melt fraction (a, left) and surface temperature (b, right). The OLR is only indirectly related to mantle melt fraction through its relationships with surface temperature and outgassed volatile inventories, which are mediated by the vertical atmospheric temperature structure and composition. The ASR for young Earth is indicated on the left panel, for a circular 1 AU orbit around a 50 Myr Sun (Baraffe et al., 2015). All models’ OLR curves show similar qualitative behaviours because thermal emission scales super-linearly with temperature. For a greybody emitter, we should expect OLR to behave as a log-linearly decreasing line in panel b, following the Stefan-Boltzmann radiation law $\varepsilon\sigma T^{4}$.

The OLR initially exceeds ${10}^{5}\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for all of these Nominal Earth simulations, except MOAI, with initial surface temperatures above 2500 K that correspond to deep magma oceans. High OLR values are achieved, even for these optically-thick atmospheres extending far above the magma ocean surface, because temperatures at the $\tau=1$ atmospheric photosphere layer are large (Figure 11). At initially high temperatures, PROTEUS finds an $\mathrm{OLR}>${10}^{5}\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$$ because its atmosphere sub-module determines that deeper, near-surface regions are convectively stable; energy is transported by radiative diffusion at a shallower lapse rate than would be supported by the corresponding adiabat in these regions (Selsis et al., 2023; Cmiel et al., 2025). Since the atmospheric temperature gradient is steeper, high photospheric temperatures and OLR values are achieved (Nicholls et al., 2025c, 2026b).

In comparison, MOAI calculates a substantially smaller initial OLR relative to the other models, when compared at a given surface temperature (Figure 8b). This lower OLR arises because MOAI sets the atmosphere $T(p)$ structure with a \chH2O pseudoadiabat and 200 K stratosphere, which results in a cooler photosphere. MOAI also predicts a highly-reduced early magma ocean which outgases \chCO-dominated atmospheres (Figure 6).

The simulations agree that OLR decreases as the planet cools and the atmosphere thickens (Figure 8b). Compositional evolution causes deviation from a simple log-linear relationship between OLR and surface temperature – especially notable for MOAI and PROTEUS (orange line in Figure 8b).

NEONGOOEY, GOOEY, LINCS, and PlanAtMO find that OLR becomes approximately constant as a function of surface temperature, for surface temperatures less than 2000 K. In these cases, OLR approximates the Simpson-Nakajima runaway greenhouse limit (dash-dot black line), except GOOEY. This OLR limit arises from their modelling of saturated upper atmosphere layers for pure-steam compositions (Nakajima et al., 1992; Ingersoll, 1969). The OLR calculated by GOOEY asymptotically tends towards $797\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ since GOOEY imposes an isothermal stratosphere, held constant at the radiative skin temperature; once the deep adiabatic region has cooled sufficiently, the $\tau_{\mathrm{LW}}=1$ photosphere layer sits within the constant-temperature stratospheric region, that determines the OLR. MOAI does not resolve the Simpson-Nakajima limit at cooler temperatures, despite also following a pseudoadiabatic temperature structure, because its radiative transfer scheme does not incorporate the effects of \chH2O collision-induced absorption (Lichtenberg et al., 2021; Zahnle et al., 1988). PACMAN maintains a comparatively high OLR, above the Simpson-Nakajima limit, due to the low atmospheric water content (Figure 3)

Our understanding of the models’ different approaches to parametrising planetary climates and radiative transfer can explain the calculated relationships between OLR and surface temperature in Figure 8b. The models’ different cooling behaviours reflect that both atmospheric composition and its temperature structure play a role in the cooling of magma ocean surfaces.

Models which exhibit comparable behaviour in OLR-$\mathrm{T_{surf}}$ space (Figure 8) are not consistent with the three compositional groups identified from Section III.3. For instance, the predictions from NEONGOOEY and PlanAtMO radiative transfer modelling (Figure 8b, yellow and pink lines) are in good agreement, correspondent with their similar predictions of outgassing behaviour (Figure 3). If we compare these models’ OLR curves at at a $\sim 95$% melt fraction (Figure 8a), the outgoing radiation flux from GOOEY is comparable to PACMAN and NEONGOOEY (red, blue, pink lines). NEONGOOEY predicts a surface temperature similar to GOOEY at this high melt fraction, yet features a \chCO2-dominated atmospheric composition, despite having comparable OLR.

On the other hand, PROTEUS shows higher OLR than both NEONGOOEY and PlanAtMO, despite similar compositional predictions, due to its radiative-convective climate calculation.

Using Figure 3, PACMAN and MOAI were classified as having similar compositional behaviours, yet PACMAN calculates the highest initial OLR and MOAI the lowest (Figure 8b, blue and purple lines). These models converge in OLR-$\mathrm{T_{surf}}$ space under cooler conditions, which follows from MOAI’s adoption of a cool 200 K isothermal stratosphere.

Clearly, models’ different approaches to climate and radiation modelling have a stronger impact on the OLR variance than the modeled evolution of atmospheric composition, because curves for substantially different compositions overlap. The models’ OLR-$\mathrm{T_{surf}}$ curves are largely intertwined, as a whole. Choices to neglect or incorporate key opacity sources (e.g. collision-induced absorption continua) lead to calculated OLR values spanning multiple orders of magnitude between the models.

The relative arrangement of OLR curves differs between Figure 8’s two panels because of the complex relationships mapping mantle melt fraction to magma ocean surface temperature. For example, PROTEUS has the lowest OLR of all models when compared at 100% melt fraction, but the second highest OLR of all models when compared at $T_{\mathrm{surf}}=3000$ K. This apparent discrepancy is physical, and is discussed in the following section.

Melt fraction and surface temperature are fundamentally different, but often correlated, physical quantities (Abe, 1997). Mapping these variables depends on their definition (Stixrude, 2014), the modelled mantle temperature profile (Bower et al., 2018), melting solidus and liquidus curves (Wolf & Bower, 2018; Labrosse et al., 2015), and mineralogical properties (Maurice et al., 2023; McDonough & Sun, 1995; Ghiorso et al., 2002). Figure 9 presents three diagnostic variables of mantle geodynamics, as functions of volumetric mantle melt fraction, from our Nominal Earth simulations: (a) surface temperature, (b) characteristic mantle viscosity, and (c) characteristic radius of solidification.

Figure 9a largely explains the differences in simulated Earth solidification timescales presented by Figure 2a. PROTEUS is systematically the ‘slowest’ model to cool, while PlanAtMO and MOAI are slightly more rapid. This line-up is reflected by surface temperature-melt fraction relationships shown in Figure 9a, where PROTEUS has the lowest $T_{\mathrm{surf}}$ at high melt fractions, followed by PlanAtMO and MOAI with slightly higher temperatures. Atmospheric temperatures and energy fluxes – which determine the rate of planetary cooling – are set by chemical composition, incoming stellar radiation, and upwelling surface emission. These quantities are directly tied to surface temperature but not directly to the mantle melt fraction. So, lower surface temperatures at high melt fractions correspond to lower outgoing fluxes (OLR) and these predictions of slower magma ocean cooling.

A parcel of silicate cooling from liquid to solid is subject to complex microphysics. Crystals initially form within a largely-molten fluid; later, melt becomes trapped within the crystalline structure of a largely-solidified mantle (Petford, 2003; Costa et al., 2009). A transition in rheological properties (e.g. viscosity) occurs between 20–60% melt fraction: firstly, a solid-liquid transition near $\Phi\approx 50\%$, secondly a connectivity transition near $\Phi\approx 7\%$ (Costa et al., 2009; Scott & Kohlstedt, 2006). All CHILI models except LINCS analytically parametrise a single rheological transition centred around a ‘critical melt fraction’ between 20 and 40%, depending on the model, whereas LINCS defines the transition radius where the local temperature equals the solidus temperature (0% local melt fraction). Figure 9b plots the radial location $R_{\mathrm{RF}}$ of the rheological transition front versus whole-mantle melt fraction $\Phi$, calculated by each model. This can be interpreted as a measure of radial melt distribution, given some bulk amount of molten mantle material. Figure 9b shows that $R_{\mathrm{RF}}$ and $\Phi$ are predicted to be strongly correlated by each model, as expected, since all scenarios exhibit bottom-up crystallisation modes and adopt qualitatively comparable melting curves (Section A).

LINCS deviates from the other models in $R_{\mathrm{RF}}-\Phi$ phase space by maintaining a deep magma ocean throughout the mantle, until the whole mantle melt fraction ($x$-axis) approaches its rheologically critical value, $\approx 50\%$— a stagnation driven chiefly by a deep adiabatic temperature profile nearly parallel to the solidus, as well as the definition of the transition radius. Its transition radius then increases abruptly to catch up with the other models.

All models present a bottom-up mode of mantle crystallisation: the rheological transition (Figure 9b) begins at the core-mantle boundary and progresses radially upwards as crystal cumulates settle (Labrosse et al., 2015; Solomatov & Stevenson, 1993). Fractional or batch (equilibrium) crystallisation behaviours depend on the cumulate sizes formed during crystallisation and the magma ocean flow regime (Labrosse et al., 2015; Bower et al., 2022). These models generally consider fractional crystallisation behaviours, except MOAI, which considers batch crystallisation. None of these simulations correspond to a basal magma ocean scenario.

Close inter-model correspondence in $R_{\mathrm{RF}}-\Phi$ space (Figure 9b) contrasts with the substantial differences in $T_{\mathrm{surf}}-\Phi$ space (Figure 9a). In the case of PROTEUS, a key difference arises from its mixing-length treatment of mantle dynamics, where upward energy transport of heat by convection is partially offset by downward transport of heat by the settling of cumulates and the release of gravitational potential energy (Bower et al., 2018; Hirschmann, 2012). Convective suppression causes the mantle temperature profile to deviate from an adiabat, yielding shallower geothermal lapse rates and lower surface temperatures. Additionally, the liquidus adopted by PROTEUS has the lowest temperature at the core-mantle boundary of all models (Figure 10a), so its simulations of Nominal Earth’s deep mantle must attain cooler temperatures than other models before solidification begins – further promoting longer solidification timescales.

Figure 9c shows that PROTEUS (orange region) calculates higher effective mantle viscosities $\eta$ compared to other models, for a given melt fraction. Compared to the other models, the orange PROTEUS curve spans a wider range of viscosity values because its interior geodynamics module (SPIDER; Bower et al., 2018) calculates 1D mantle viscosity profiles, which depend on the melt fraction of each layer, while other models solve for a scalar value representative of the whole mantle (Costa et al., 2009; Bottinga & Weill, 1972). However, even for the initially fully-molten magma ocean conditions, e.g. where $\Phi\approx 95\%$, the models predict a range of viscosity values which vary across orders of magnitude (panel b). This reflects their different viscosity parametrisations. These model differences can be readily understood by considering the boundary layer parametrisation of Rayleigh-Bernard mantle convection (Solomatov, 2015; Schaefer et al., 2016; Lebrun et al., 2013). In boundary layer theory, the convective heat flux scales as a function of the Rayleigh number,

The Rayleigh number, Ra, characterises the relative importance of buoyancy-driven convective motion compared to the resistive effect of thermal and momentum diffusion in a fluid (Turcotte & Schubert, 2002). Ra quantifies the mantle’s transition from an initially turbulent regime towards a laminar flow regime. The empirical exponent $n$ is often taken to be $2/7$ in the hard turbulence and $1/3$ in the soft turbulence regimes (Turner & Campbell, 1986; Turcotte & Schubert, 2002). The CHILI models’ parametrisations of mantle geodynamics (Section II.2) generally use $n=1/3$ and a scaling for Ra which is inversely proportional to mantle viscosity:

where $T_{\mathrm{pot}}$ is the mantle potential temperature and $\eta$ is the melt fraction or temperature-dependent dynamic viscosity of the mantle. For example, GOOEY adopts an Arrhenius-type viscosity function $\eta\propto\exp[-E_{a}/(R_{\mathrm{gas}}T_{\mathrm{pot}})]$ (Schaefer et al., 2016). PROTEUS uses an Arrhenius-type function with a depth-dependent activation volume (Bower et al., 2018). Thus, higher viscosities make mantle convection less efficient and inhibit heat transport out of the deep interior (Bower et al., 2022; Turner & Campbell, 1986; Solomatov, 2015). That PROTEUS calculates higher viscosities is correspondent with it finding longer cooling timescales in Figure 2.

Figure 9c shows that our models predict only late-stage increases in viscosity to values above 1 Pa s, at melt fractions $<5\%$. These differing behaviours are partially geometric: the volume-averaged melt fraction ($x$-axis) is weighted towards molten upper mantle layers which have larger volumes. However, it is clear from cooling timescales in Figure 2 and geodynamical diagnostics in Figure 9 that different approaches to parametrising the microphysical processes of silicate crystallisation correlate with predicted cooling timescales differing by orders of magnitude.

Unmodelled physics, such as the rheological effects of volatiles dissolved within the magma ocean fluid, would suggest that these calculated magma ocean viscosities are upper-limits on real behaviours. Inter-model agreement of viscosities far below that of a purely solid phase (dash-dot line) across a wide range of melt fractions is therefore a robust prediction from these evolutionary simulations. Real planetary magma oceans are expected to be highly inviscid.

All of these models each incorporate a type of parametrised, thermally conductive boundary layer (CBL) at the atmosphere-mantle interface. This CBL is different to the tectonic crust of present day Earth, and is intended to represent a potential floatation ‘skin’ atop the magma ocean induced by radiative cooling (Elkins-Tanton, 2008; Abe & Matsui, 1985). The CBL cannot be spatially resolved by these models but can, in principle, inhibit cooling by decreasing the calculated surface temperature (and thus net thermal emission) relative to the mantle potential temperature (Bower et al., 2018; Spaargaren et al., 2020). GOOEY, NEONGOOEY, and PACMAN predict substantial boundary layer thickening at later stages of Nominal Earth simulations (Appendix C). CBL growth to $\sim 1$km thickness occurs during the final stages of NEONGOOEY and GOOEY simulations, when mantle melt fractions are less than 0.8%, thereby having only a small opportunity for inhibiting magma ocean cooling. Yet, PROTEUS simulations impose a constant small CBL thickness of $1$cm while predicting comparable Nominal Earth solidification timescales to GOOEY (Figure 1). We conclude that the presence of a conductive boundary layer atop primordial magma oceans would have a negligible effect on their thermal evolution, if it were to form (Elkins-Tanton, 2008), although it could become relevant after a transition to post-magma ocean tectonic regimes (Korenaga, 2013; Tackley, 2000; Meier et al., 2026; Maurice et al., 2017; Gillmann & Tackley, 2014).

All models are able to reproduce Earth and Venus solidification time-scales consistent with empirical constraints (Figure 1). The simulated termination ages for an early magma ocean in our Nominal Earth scenario occur long before lunar formation, at ages between $50\sim 200$ Myr (Borg et al., 2014; Barboni et al., 2017), allowing for multiple distinct magma ocean epochs during the Hadean. This is consistent with indications from measurements of highly fractionated ³He/²²Ne isotope ratios from mid-ocean ridge degassing (Tucker & Mukhopadhyay, 2014). Additionally, the solidification time-scales for all Earth cases (Nominal case in Figure 1, and at various H+C inventories in Figure 2) respect the radiometrically estimated ages of Hadean mafic intrusions dated to have solidified within 300 Myr of Earth’s formation (Dalrymple, 2001; Sole et al., 2025)

We also expect that Venus would take longer to solidify than Earth, all else equal, given its higher insolation (Hamano et al., 2013; Krissansen-Totton et al., 2021b). Our models reproduce this expectation (Figure 1).

However, there are quantitative differences between the models’ simulated evolution of Nominal Earth and Venus. The evolutionary timescales for reaching sequential melt-fraction thresholds (Figure 2) vary across orders of magnitude, reflecting our different modelling assumptions and approximations. We have identified the key assumptions which shape these different timescales.

Firstly, atmospheric composition strongly regulates planetary cooling through gas radiative opacities blanketing the interior, inducing a greenhouse effect (Section III.3). The models’ approximations of atmospheric chemistry are thereby reflected in the cooling timescales through climate modelling. In particular, models’ approaches to calculating vertical atmospheric temperature structure determines the planet’s net outgoing radiation flux (Figure 8). However, our models agree that a planet’s post-formation hydrogen budget strongly dictates magma ocean cooling rates (Hamano et al., 2015; Nicholls et al., 2024).

Secondly, since these modelled atmospheric compositions are directly calculated from mantle oxidation state, interior redox processes strongly shape atmospheric chemistry. Outgassed \chCO2/\chCO ratios vary as a function of degassing temperature and mantle ferric-ferrous iron partitioning (Section III.3), so fully-coupled treatments of planetary interior and atmospheric evolution are critical to self-consistent evolution modelling. Atmospheric escape remains a substantial point of model disagreement (Figure 7), although our Venus simulations indicate that volatile element fractionation plays a minor role in regulating magma ocean thermal evolution (Section III.4).

Thirdly, mantle temperature structure, energy transport processes, and assumed phase-state boundaries (i.e. melting curves) substantially influence when and where planetary interiors are (semi-)molten. Model inter-comparison is best undertaken at states of similar melt fraction, rather than surface temperatures, because model $\Phi-T_{\mathrm{surf}}$ relations differ substantially. In addition to our models’ adoption of different melting curve parametrisations (Appendix A), mantle solidus and liquidus boundaries are also sensitive to multiple parameters and physical effects not considered by these codes (Katz et al., 2003; Ghiorso et al., 2002; Baumeister et al., 2025). For example, the incorporation of dissolved water (as OH^- ions) into the molten silicate is effective at disrupting crystal polymerisation (Katz et al., 2003; Whittington et al., 2000). These dissolved volatiles have the effect of ‘solidus depression’, which readily reduces the solidus temperatures by $\sim 100$s of kelvin, depending on the volatile content, to further sustain molten states (Elkins-Tanton, 2008). Figures 2 and 9 show that differences in these existing thermodynamic modelling choices lead to cooling timescales across multiple orders of magnitude, so careful future consideration of additional physical sensitivities in mantle melting properties is advised (Thompson, 2026).

That our models agree upon bottom-up solidification scenarios without basal magma oceans could be viewed as physical robustness. However, it largely arises from all models’ parameters being tuned to reproducing a bottom-up crystallisation mode by their choice of melting curves (Elkins-Tanton, 2008; Labrosse et al., 2015; Solomatov & Stevenson, 1993). More substantial inter-model deviations would be expected had they adopted different solidii, or considered a diverse range of crystal cumulate sizes that have various proclivities for settling-out of turbulent magma ocean flow (Spaargaren et al., 2020; Bower et al., 2018; Lichtenberg et al., 2023).

The lifetime histories of Earth and Venus raise numerous open questions. For example, there are substantial uncertainties in Earth’s current volatile budget (Peslier et al., 2017; Lichtenberg et al., 2019; Marty & Genda, 2025) and the abundance of light element in its metallic core (Labrosse et al., 1997; Schoenberg et al., 2002; Shahar et al., 2026). Synergistic application of planetary formation models and coupled interior-atmosphere simulations can, in principle, constrain the elemental compositions of these reservoirs by tracking volatile budgets. This desire also holds for Venus (Hamano et al., 2025), Mars (Drilleau et al., 2026; Maurice et al., 2017), and exoplanet scenarios (Nicholls et al., 2026a; Barth et al., 2021). However, our model intercomparison (Section III) suggests that physical parameters to be inferred (e.g. a planet’s initial H content) are largely overprinted by the substantial differences in predicted planetary evolution behaviours which arise from our modelling choices and parametrisations.

We have shown that simulated outcomes depend on multiple physical and chemical processes, which exhibit different levels of sensitivity in their relative importance and modelled parametrisations, and mutually interfere in either constructive or detractive modes. Our controlled intercomparison of planetary evolution models has benchmarked these processes in the context of Venus and Earth. It is critical to identify the physical and chemical processes which (i) must be included in models because they exert leading-order control over planetary composition and evolution, and (ii) must be incorporated into models accurately. Ideally, by focussing on these processes, we can work towards a point where evolutionary models agree on planetary cooling timescales and compositions. Table 2 presents the key physical parametrisations adopted by our models and, on the basis of this intercomparison, categorises these processes’ relative importance in controlling planetary thermo-compositional evolution.

We can better calibrate our models by adopting additional measurement constraints. In this paper, we have tentatively used estimates for CAI ages (Baker et al., 2005), Earth’s age at core-differentiation (Dalrymple, 2001), and the ages of the oldest mineral samples (Sole et al., 2025; Mojzsis et al., 2001), to bracket the valid range of timescales on Earth’s early thermal evolution (Figures 1 and 2). All models here respect these age constraints, but these leave orders of magnitude flexibility in cooling timescale, and do not directly constrain compositions or climate states.

The simulation of other governing physics within these codes may enhance model accuracy (Lichtenberg et al., 2025; Nicholls, 2026b; Schaefer & Elkins-Tanton, 2018). Meanwhile, further model developments to predict additional historical fingerprints and proxies would be valuable as points for calibration against reality.

For example, future iterations of these codes could calculate D/H isotopic fractionation due to hydrodynamic escape; comparison of models’ predicted D/H ratios against each other (Schaefer et al., 2024), and against empirical measurements from Earth (Peslier et al., 2017) and Venus (Grinspoon, 1993), would inform us about the accuracy of their escape implementations. Another opportunity for model intercomparison and empirical calibration is through parametrisation of gas-melt and melt-crystal noble gas partitioning into these codes. Earth and Venus measurements of argon (O’Rourke & Korenaga, 2015; Kaula, 1999) and xenon (Rzeplinski et al., 2022) fractionation are suggestive of particular magma ocean crystallisation scenarios. Yet, the CHILI models do not homogeneously quantify metrics for comparison against these measurements (Krissansen-Totton et al., 2021b; Gilbert-Janizek et al., 2026).

No model is ever ‘correct’, but we may strive towards sufficient model robustness and complexity in order to understand and explain measurements of the natural world (Box, 1976; Sohl et al., 2024). Here, we discuss the potential for additional, unmodelled processes to shape planetary compositions and structures. These processes could be considered for incorporation into the CHILI codes and similar modelling frameworks at future stages of their development. Doing so may resolve the disagreements presented here (Section IV.1) and enable their valid application to other physical regimes.

Figure 9 highlights that these models’ different approaches to simulating mantle temperature structures and mantle rheological properties can partially explain their different predicted magma ocean cooling timescales. Ultimately, these differences in model approach reflect the complex solidification microphysics which regulates crystallisation, cumulate entrainment, and Stokes settling (Wordsworth et al., 2018; Lichtenberg, 2021; Schaefer et al., 2024). Future model developments should carefully incorporate robust parametrisations of these effects, which are, in some scenarios, also expected to generate basal magma ocean outcomes (Boukaré et al., 2025, 2018; Labrosse et al., 2015; Lark et al., 2026). This potential outcome is an important consideration, since we have shown that deep volatile storage regulates outgassed atmospheric compositions (Section III.3).

We have identified important differences between CHILI models’ predicted Earth and Venus evolution pathways. Although models will necessarily adopt simplifications of reality, given finite computational resources, their parametrisations reflect our epistemological uncertainties the underlying physics. For example, various atmospheric compositions are generated due to different treatments of solid-phase melt trapping (Section III.3). Yet, the partitioning coefficients of volatiles into mineral phases remain poorly constrained by experiments (Schaefer et al., 2024). Differences in models’ predicted OLR motivate the adoption of accurate radiative transfer parameterisation, which are especially limited by unconstrained \chH2O continuum absorption coefficients (Abel & Frommhold, 2013; Pierrehumbert, 2010).

All of these codes treat the planet’s metallic core as compositionally inert and leave it spatially unresolved. Yet, core fluid dynamics, cooling, magnetism, and the core-mantle partitioning of volatiles are rich fields of study (Labrosse et al., 1997; Lay et al., 2008; Stevenson, 2001). The deep core domain of planetary interiors could be incorporated into these interior-atmosphere evolution models (Huang & Dorn, 2025; Garcia et al., 2026). This could, for example, permit model predictions regarding the impact of dynamo-generated magnetism on atmospheric escape (Gülcher et al., 2020). Magnetic effects are suggested to enhance mass loss from Venus at the present day (Lammer et al., 2020; Zahnle & Kasting, 1986).

The CHILI protocol specifies fixed orbital parameters to limit intercomparison scope. However, planetary and satellite orbital dynamics are strongly coupled planetary interiors (Barnes, 2017; Korenaga, 2025; Farhat et al., 2022). Lunar tides and resonant states offer one explanation for the in-situ evidence of lunar surface remelting (Nimmo et al., 2024; Chen & Nimmo, 2016; Farhat et al., 2025). Theoretical models can resolve the important coupling between tidal heating, orbital evolution, and thermal regulation by overlying atmospheres (Zahnle et al., 2015; van Dijk et al., 2026). Theoretical modelling of tidal heating within multi-planetary exoplanet systems also suggests that magma oceans may be sustained under temperate irradiation regimes (Hay & Matsuyama, 2019; Herath et al., 2024; Unterborn et al., 2018; Pu & Lai, 2019). Incorporation of tidal heating effects and orbital migration into planetary evolution models would be valuable, but may necessitate a careful intercomparison (Sabadini et al., 2016).

The CHILI protocol also specifies clear-sky climate modelling. However, clouds strongly determine Earth’s present day climate and weather (Pierrehumbert, 2010, 2002). Venus’ photochemical sulfur clouds are the source of its high $\approx 75\%$ Bond albedo (Titov et al., 2018), and are crucial to suggested greenhouse feedback effects early in its lifetime (Krissansen-Totton et al., 2021b; Hamano et al., 2025). The presence of clouds may have determined whether Venus once hosted surface water (Turbet et al., 2021; Kitzmann et al., 2010; Pollack, 1971). More generally, aerosols can induce a greenhouse or anti-greenhouse effect depending on particle size distributions, particle shapes, compositions, and their vertical distribution (M. S. Marley & A. S. Ackerman, 2014; Helling, 2019; Lodge et al., 2026). Aerosols could extend or foreshorten the magma ocean solidification timescales predicted by the CHILI models here (Piette et al., 2023). The ‘problem’ of aerosols can be approached at various levels of modelling complexity (Helling, 2019), so it remains uncertain how best to study their role in deep planetary histories and resultant impact on present climates. MOAI incorporates Rayleigh scattering but neglects key collisional continua (e.g. \chH2O-\chH2O; Lichtenberg et al. (2021)), and predicts Earth and Venus magma ocean solidification timescales intermediate to the other CHILI models (Figures 1 and 2). The CHILI models also adopt different approaches to calculating atmospheric temperature structure, which, alongside atmospheric composition, sets the planetary cooling rate and potential for cloud formation (Cmiel et al., 2025; Nicholls et al., 2025c; Peng & Valencia, 2024; Janssen et al., 2026). Overall, understanding the behaviours emergent from the intimate competition between thermodynamic and radiation processes, which we have shown to play a role in setting magma ocean cooling rates (Section III.5), requires a dedicated intercomparison of exoplanet climate models. The ongoing COD-ACCRA and MALBEC intercomparison projects, also part of CUISINES, will directly address these uncertainties (Villanueva et al., 2024; Chaverot et al., In prep.).

The discovery of sub-Neptune and super-Earth exoplanet populations overturned the picture of planetary formation and evolution established by the Solar System planets (Bean et al., 2021; Valencia et al., 2006). Yet, current telescopes (e.g. JWST) and upcoming observational surveys (e.g. PLATO) can only probe exoplanets’ bulk compositions and upper-atmosphere chemistry. Understanding the true nature of these new planetary populations requires comprehensive predictive and interpretive modelling approaches (Lichtenberg et al., 2025). However, planetary evolution models have largely targeted Earth-sized planets, which remain the focus for CHILI parts I, II, and III. Modelling of sub-Neptune magma oceans and thermal evolution may require a careful handling of additional physics. For example, sub-Neptunes’ potential for interior-atmosphere (im)miscibility and high temperature supercritical fluid layers that necessitate the extension of thermodynamic data to extreme pressures and temperatures (Rogers et al., 2025; Young et al., 2025; Wolf & Bower, 2018; Attia et al., 2026). The upper-atmosphere compositions probed by telescopes are also subject to disequilibrium chemical effects, potentially requiring the consideration of complex upper-atmospheric photochemistry within these coupled interior-atmosphere modelling frameworks (Nicholls et al., 2023; Tsai et al., 2021; Yu et al., 2021).

Ultra-short period exoplanets are highly observable with JWST through secondary eclipse and phase curve measurements (Hammond et al., 2025; Teske et al., 2025; Demory et al., 2016). The magma oceans on these ‘lava worlds’ may be permanently sustained by intense stellar irradiation (Lichtenberg & Miguel, 2025), so application of planetary evolution and climate modelling to these regimes can enable contemporary insight into the primordial environments of all rocky worlds (Janssen et al., 2026; Nguyen et al., 2024). CHILI Part II and Part III will address the opportunity afforded by lava worlds, to test the validity of coupled interior-atmosphere modelling under intense irradiation and identify the key physics shaping their observations. These planets’ exposure to intense irradiation may lead to substantial atmospheric escape, in excess of that found for our Earth and Venus scenarios, and thereby modulate their thermal evolution.