{'fraction_non_alphanumeric': 0.03377157234487402, 'fraction_numerical': 0.023707431163400544, 'mean_word_length': 4.675113122171946, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low volcanic outgassing rates for a stagnant lid Archean Earth with graphite-saturated magmas

Volcanic gases supplied a large part of Earth's early atmosphere, but constraints on their flux are scarce. Here we model how C-O-H outgassing could have evolved through the late Hadean and early Archean, under the conditions that global plate tectonics had not yet initiated, all outgassing was subaerial, and graphite was the stable carbon phase in the melt source regions. The model fully couples numerical mantle convection, partitioning of volatiles into the melt, and chemical speciation in the gas phase. The mantle oxidation state (which may not have reached late Archean values in the Hadean) is the dominant control on individual species' outgassing rates because it affects both the carbon content of basaltic magmas and the speciation of degassed volatiles. Volcanic gas from mantles more reduced than the iron-w\"ustite mineral redox buffer would contain virtually no CO2 because (i) carbonate ions dissolve in magmas only in very limited amounts, and (ii) almost all degassed carbon takes the form of CO instead of CO2. For oxidised mantles near the quartz-fayalite-magnetite buffer, we predict median CO2 outgassing rates of less than approximately 5 Tmol/yr, still lower than the outgassing rates used in many Archean climate studies. Relatively weak outgassing is due in part to the redox-limited CO2 contents of graphite-saturated melts, and also to a stagnant lid regime's inefficient replenishment of upper mantle volatiles. Our results point to certain chemical and geodynamic prerequisites for sustaining a clement climate with a volcanic greenhouse under the Faint Young Sun.


Introduction
The Hadean-Archean atmosphere cloaked the only environment in the universe known to support the origin of life. This atmosphere would have been fed by volatile chemicals outgassed from the interior (Holland, 1984;Kasting et al., 1993;, as well as extraterrestrial impacts  and perhaps a lingering, degassed magma ocean atmosphere ). Yet hardly any observable record of the Archean atmosphere remains. Our knowledge surrounding its makeup rests on modelling. Mantle outgassing models lend boundary conditions to the processes controlling the atmospheric composition and planetary climate. However, predictions from these models are very much entwined with certain assumptions-in particular, the tectonic regime and oxidation state of the mantle. Archean outgassing scenarios hinging on mobile plates and an oxidised mantle form the basis of contemporaneous climate studies (e.g., Sleep and Zahnle, 2001). This work offers an alternative to these recurring premises: a fully-coupled numerical estimation of the volcanic outgassing rates of CO 2 , CO, H 2 O, and H 2 , for a stagnant lid Earth under a range of mantle oxidation states.
The oxidation state of the upper mantle controls the chemical speciation of volcanic gas. Increasing the mantle's oxygen fugacity, f O2 , raises the ratio of oxidised to reduced volatile species in the melt phase: [CO 2 (Holland, 1984). As the melt rises adiabatically through conduits, its f O2 re-equilibriates with the decreasing pressure; gas abundances at the surface are also linked to the mantle source region (Kasting et al., 1993). On modern Earth, for example, the dominant species in volcanic gas are H 2 O and CO 2 , in line with a mantle source f O2 in equilibrium with the quartz-fayalitemagnetite (QFM) mineral buffer (Holland, 1984).
Besides the effect of f O2 on volatile speciation, there can be a compounding redox control on the partitioning of carbon into the melt. This effect becomes relevant when graphite is the stable carbon phase in the upper mantle. During melting, carbonate ions are produced from a reaction of graphite with oxygen via CO 2 (Holloway et al., 1992). The total carbon outgassing flux depends on how much carbon (as CO 2-3 ) can dissolve into the melt.
Although these two redox effects are well-known, the value itself of the Archean mantle f O2 remains weakly constrained at best. Immediately after core formation, the magma ocean must have been reduced (Wood et al., 2006). Meanwhile, scant geological evidence suggests that by ∼3.8 Ga at the latest, the uppermost mantle was already about as oxidised as today (Canil, 1997;Delano, 2001;Trail et al., 2011;Nicklas et al., 2019;Armstrong et al., 2019). Yet uncertainty persists on the timing of this oxidation (see section 2.5.3). Therefore, to cover the possible spread of redox states, we investigate outgassing fluxes under a wide range of mantle f O2 , especially in proportion to other unknown parameters such as the mantle volatile content and thermal state.
Whilst this study is set up to explore many redox scenarios, it will be limited to one tectonic regime. For this we consider the end-member case of stagnant lid convection (Debaille et al., 2013). The possibility of a stagnant lid regime appears underexplored in early Earth climate studies, despite the ambiguity around plate tectonics' operation through the Archean (e.g., Brown et al., 2020); as such we hope to complement other work assuming a mobile lid. This modelling decision is important because mobile lids may be associated with much higher outgassing rates than stagnant lids (Noack et al., 2014(Noack et al., , 2017. The mechanistic explanations are not straightforward-being a combination of competing effects discussed later in section 4.1.5-and would benefit from dedicated work in the future. By numerically coupling mantle convection, redox-dependent melt partitioning of volatiles, and redoxdependent speciation, and by investigating a broader range of mantle oxidation states, we can advance the current literature dedicated to C-O-H outgassing from solid-state mantles of rocky planets in the solar system and elsewhere (O'Neill et al., 2007;Grott et al., 2011;Noack et al., 2014Noack et al., , 2017Tosi et al., 2017;Dorn et al., 2018;. Previous studies have tended to assume an oxidising mantle, akin to that of present-day Earth-exceptions include Grott et al. (2011) for the reduced early Mars, and Tosi et al. (2017), , and Liggins et al. (2020) for fiducial rocky exoplanets. Further, numerical convection carries the advantage of resolving the local temperatures and pressures of melt parcels, to which f O2 is very sensitive.
An approach parallel to the one taken here is to parameterise outgassing fluxes in terms of melt production (e.g., Sleep and Zahnle, 2001;Papuc and Davies, 2008;Kite et al., 2009;Kadoya and Tajika, 2015;Foley and Driscoll, 2016;Charnay et al., 2017;Krissansen-Totton et al., 2018;Foley and Smye, 2018;Krissansen-Totton and Catling, 2020). In parameterised models, mantle melting rates are often scaled to an estimate of the interior heat flux. Past outgassing rates can then be calibrated to present-day estimates. These types of outgassing parameterisations may imply higher outgassing fluxes than do numerical convection models for two reasons. Firstly, heat flow scalings extrapolated backwards in time will necessarily lead to higher-than-modern outgassing rates because parameterised thermal history models predict high heat flow early in a planet's evolution (e.g., Turcotte, 1980). Secondly, 1D models tend to overestimate melt fractions early in a planet's thermal evolution, since they cannot spatially resolve depletion upon partial meltinglocal depletion would inhibit further local melting for the same temperature. In addition to these systematic effects, applying a single scaling relationship assumes no change in mantle oxidation state or tectonic regime.
Planetary outgassing rates require attention because they lie at the core of at least two unsolved Archean mysteries. The first is that reduced gases such as H 2 and CO may be necessary for prebiotic chemistry. The production of HCN requires a reducing atmosphere with C/O ≥ 1, for example (Rimmer et al., 2019). How-ever, whether H 2 and CO could have comprised a significant part of volcanic gas at the time remains debated (e.g., . The second mystery is that sufficiently large partial pressures of greenhouse gases seem required to incubate early Earth under lower solar luminosity (Charnay et al., 2020) and make sense of liquid surface water (Sagan and Mullen, 1972). Although no solution to this paradox has emerged as the likely favourite, most atmospheric evolution scenarios proposed so far rely on CO 2 outgassing fluxes at least as high as the present day (e.g., Sleep and Zahnle, 2001;Wordsworth and Pierrehumbert, 2013;Kanzaki and Murakami, 2018). Neither of these two driving questions may ever be resolvable without definite constraints on the timing of mantle oxidation, or plate tectonics' initiation. Nevertheless, this study, by more rigorously linking stagnant lid outgassing fluxes to the unknown oxidation state, might help contextualise these problems within Earth system evolution.

Methods
We trace the movement of volatiles through the mantle to the atmosphere for a stagnant-lid Earth, from just after the moon-forming impact to around the first rock record. Our model is similar to the one described in , who additionally calculated gas solubility in the melt (relevant for the higher outgassing pressures expected on exoplanets greater than an Earth mass), but coupled partitioning and volatile speciation a posteriori within the convection simulations, rather than directly.
To simulate mantle dynamics, we use a 2D convection model, which is detailed in Noack et al. (2016). The heavily-benchmarked model couples solid-state convection in the mantle with spatially-resolved melt fractions and volatile partitioning in the melt (Holloway et al., 1992;Katz et al., 2003;Grott et al., 2011). This is further coupled with a chemical speciation model (French, 1966;Holloway, 1981;Fegley, 2013;Schaefer and Fegley, 2017). At each time step the model produces outgassed masses of H 2 O, H 2 , CO 2 , and CO. We simulate 700 Myr of convection, with time-zero nominally corresponding to 4.5 Ga.

Numerical convection model
The convection code solves the conservation equations of mass, momentum and energy in the rocky mantle by using a particle-in-cell method. We use the 2D spherical annulus method in Hernlund and Tackley (2008) to divide the mantle into a mesh of cells with an average lateral and radial resolution of 50 km between grid points. The radial resolution coarsens linearly with depth, varying from 10 km at the surface to 90 km at the core-mantle boundary. To reduce computational costs, we model one quarter of the 2D annulus, corresponding to 307 cells in the lateral direction. We found that minor changes in resolution did not affect degassing.
We model compressible convection based on the truncated anelastic liquid approximation (TALA), using pre-calculated depth-dependent profiles for thermodynamic parameters such as the mineral-dependent density, the thermal expansion coefficient and the heat capacity for an adiabatic temperature profile. We neglect any possible initiation of plate tectonics or other surface recycling mechanisms and instead allow for a stagnant lid to form.
We take into account time-dependent radiogenic heating and heat loss from the core. The radiogenic heat sources decline over time and would reach present-day Earth values of 79.5 ppb Th, 240 ppm K and 20.3 ppb U after 4.5 Gyr of thermal evolution, following Schubert et al. (2001). We did not include the partitioning of heat sources in the crust, and instead applied homogeneous radiogenic heat sources over the entire domain.
We assume a hydrated pyrolitic mantle with an Arrhenius viscosity, η, that depends on the local temperature T and pressure p, following Karato and Wu (1993): with, for the upper mantle, a prefactor A = 3.7 × 10 −19 m −r Pa −1 s −1 , grain size d = 1 mm, grain size exponent r = 2.5, activation energy E = 240 kJ mol −1 , and activation volume V = 5.0 cm 3 mol −1 , where R is the gas constant. The lower mantle rheology behaves as in Noack et al. (2017), with a viscosity prefactor and pressure-dependent activation enthalpy following Tackley et al. (2013); that is, V decreases with pressure in the lower mantle. Whereas Tackley et al. (2013) increased their viscosities by a factor of 100 to allow for numerically-feasible simulations, we do not also adopt this scale factor, so the lower mantle is less viscous compared to their study. Whilst we model a wet mantle, we do not consider local viscosity variations due to heterogeneous water contents. We do not implement an aesthenosphere, and may underestimate the convective strength locally in the upper mantle. Our initial Rayleigh number, estimated for a global average viscosity of about 10 21 Pa s, is around 5 × 10 7 . However, depending on the initial temperature profile and its local variations, the viscosity can be much lower than this-especially in the upper mantle-permitting more vigorous convection in some regions.

Melting model
The petrological model simulates partial melting of the silicate mantle. Melt forms if the local temperature is greater than the mantle solidus temperature, T solidus , at that temperature and pressure. The local melt fraction, F , is The solidus and liquidus parameterisations, given in Noack et al. (2014), are fits to experimental data for peridotite melting valid up to 15 GPa (De Smet, 1999). The liquidus temperature, T liquidus , is constant for a given pressure. The solidus temperature decreases with water content in the melt (Katz et al., 2003), and increases upon depletion due to prior melting. Melting also consumes latent heat, feeding back into the local mantle temperatures. We assume that the upper mantle melt is buoyant only at pressures of less than 12 GPa (Ohtani et al., 1995), hence we neglect any contribution from melting at greater depths. The local melt fractions add up to give the total chemical depletion of the mantle. Here, depletion is adapted to the fact that dehydration increases the solidus temperature. We recalculate the depletion at each time step using the updated solidus temperature. Upon 30% depletion of the primordial peridotite mantle, the residual mantle therefore always resembles the harzburgite composition, independently of if it having started wet or dry. Depletion is limited to 30%. We do not consider density variations in the mantle residue upon melting; that is, there is no buoyancy effect due to depletion.

Volatile partitioning model
At each time step, if the local melt fraction is greater than zero, we calculate the amounts of C and H 2 O that enter the melt phase, and subtract them from the local chemical inventories. These inventories (as well as other local information such as the depletion) are traced via mass-independent particles that flow along convective stream lines. We ensure that each cell, which is fixed in space, has between 3 and 8 particles. No local partitioning occurs once a particle's volatile inventory is emptied; if a particle is dehydrated, it remains so. As the mantle convects, however, cells can be replenished by fresh particles rising up from the lower mantle, which replace the dehydrated particles.
Carbon dissolves into the melt as CO 2-3 ; its mole fraction, X melt , depends strongly on f O2 : where K I and K II are the equilibrium constants governing, respectively, the formation of CO 2 from graphite and CO 2-3 from CO 2 , Heat capacity 1793 J kg −1 K −1 Lesher and Spera (2015) where a = 40.07639, b = 2.53932 × 10 −2 , c = 5.27096 × 10 −6 , d = 0.0267, e = −6.24763, f = 282.56, g = 0.119242, T is in K, and p is in bar (Holloway et al., 1992;Grott et al., 2011;Hirschmann and Withers, 2008). The amount of CO 2 in weight percent, here denoted χ to distinguish from mole percent, is calculated from X melt CO 2-3 , assuming thermodynamic equilibrium: where M CO 2 is the molar mass of carbon dioxide in g mol −1 and FWM is the formula weight of the melt normalised to one oxygen, equal to 36.594 g mol −1 for an example tholeiitic basaltic (Holloway et al., 1992;Grott et al., 2011). The formula weight typically varies only by few percent for different basalt compositions. Note that although (3-6) do not directly involve the mantle source concentration, χ melt CO 2 is capped at the source concentration locally.
We assume that all solid carbon exists as graphite, since we only consider shallow melting below 12 GPa. This assumption may not be valid for the most oxidised cases (IW + 4), as carbonate may become the stable carbon species for these pressures and f O2 (Stagno, 2019). However, in such oxidised conditions we already know that CO 2 is outgassed, and neglect any carbonate stability.
For water partitioning, we use the batch melting formula: where χ rock H 2 O is the local weight fraction of water in the solid mantle, F is the local average melt fraction, and D H 2 O is a partition coefficient, which we set at 0.01 based on the Ce partitioning coefficient (Michael, 1995;Katz et al., 2003).
Whilst the melt rises adiabatically towards the surface, we only allow a fixed fraction of the melt, χ extr , to reach the surface as extrusive volcanism. Only this extrusive melt contributes to outgassing into the atmosphere. The mantle volatile concentrations are small enough to not saturate (Katz et al., 2003). We do not consider submarine degassing, so virtually all gas molecules enter the atmosphere at these lower surface pressures. Our model therefore assumes complete outgassing of the melt that reaches the surface, with no residual volatiles remaining in the melt. Finally, we do not consider any possible re-equilibration or contamination as material rises through the lithosphere, and neglect whether the melt interacts with the surrounding rock or with fluids at higher pressures.

Volatile speciation model
We retrieve the masses of outgassed volatiles using the "Equilibrium constants and mass balance method" first appearing in French (1966), and used widely in other studies (Holland, 1984;Fegley, 2013;Schaefer and Fegley, 2017; assuming that the volcanic gases are in equilibrium with their magmas. These chemical reactions have respective equilibrium constants K III and K IV , which are calculated as follows: where R is the universal gas constant (8.314 J K −1 mol −1 ), X i is the mole fraction of species i, and ∆ r G 0 and ∆ r G 0 (9) are the Gibbs free energies of reaction in J mol −1 for (8) and (9) respectively-functions of temperature given in , based on Fegley (2013). Equations (10) and (11) assume surface pressure p surf = 1 bar, whilst T is taken to be the temperature the melt would have once it has risen adiabatically to the surface: where C p,melt is the heat capacity of the melt in J kg −1 K −1 , α melt is the thermal expansion coefficient in K −1 , ρ melt is the density in kg m −3 , and T melt and p melt are the local temperature in K and pressure in Pa of the melt source. Values for the thermodynamic parameters are listed in Table 1. From the model outlined above, we can calculate the outgassed mass M i in kg of each volatile species i at a given time, accumulated over all previous time steps and over all cells in the domain: where χ extr is the percentage of volcanism that is extrusive, F j,k, is the local melt fraction in a cell, V k is the volume of that cell in m 3 , m i is molar mass of species i in kg mol −1 , and ρ mtl is the local density of the solid mantle in kg m −3 . Note that χ has units of weight fraction and X has units of mole fraction. The melt volumes correspond to a 2D annulus, later scaled up to 3D by a dimensionalisation factor. The partial pressure in Pa outgassed to the atmosphere for species i is: where M tot is the sum over species of the outgassed masses from equation (13) in kg, X i,tot is the species' mole fraction, g is the acceleration due to gravity at the surface (9.8 ms −2 ), and R E is the radius of Earth (6371 km). We assume that the gas reaches the atmosphere instantaneously.

Initial conditions and free parameters
The main study involves >700 cases with random values for eight free input parameters (Table 2), identified as the main unknown parameters potentially affecting the early Earth evolution. Input values for each case are drawn from uniform distributions, with one exception. This section outlines the choice of bounds.

Choice of initial temperatures
The initial mantle temperature profile is set by four unknown parameters: planet surface temperature, T surf , initial stagnant lid thickness, D lid , initial temperature at the base of the lithosphere, T ini , and core temperature jump, ∆T c . Temperatures increase linearly from T surf to T ini , and then follow an initially adiabatic profile from below the upper thermal boundary layer to the core-mantle boundary. The temperature contrast between the base of the mantle and the core is fixed at ∆T c . We enforce an initial temperature cut-off at the solidus.
Although evidence of surface temperatures covering our age of interest is scant, oxygen isotope ratios from the 4.4 Ga Jack Hills zircons imply a surface temperate enough to host liquid oceans (Valley et al., 2002). Thus we allow T surf to vary from 0 to 60 • C. Our constant T surf is an approximation; in reality this value varies spatially and is linked to outgassing via greenhouse warming.
The temperatures in Earth's mantle after the last magma ocean stage are unknown. We adopt a range of 1750-2000 K in T ini . This study's effective distribution of T ini is not uniform because higher values of T ini were more likely to result in numerical instabilities. Values of T ini > 2000 K always led to numerical errors. However, temperatures do increase beyond this initial value due to radiogenic heating. Archean komatiite records suggest an upper mantle potential temperature of at most ∼1970 K (Herzberg et al., 2010), corresponding to upper mantle temperatures comparable with those investigated here (our T ini values refer to the actual temperatures below an initial upper thermal boundary layer, hence the potential temperatures would be slightly smaller).
The initial adiabatic profile implies an already-convecting mantle. However, at the onset of our simulations, the mantle is gravitationally stable. Hence we slightly perturb the temperature field to trigger convection. Initial temperatures are close to the solidus, so melting and outgassing begin quickly. Regardless, this temperature profile quickly starts evolving, and its initial shape is not significant to the results. The stagnant lid thickness also adjusts according to the temperature profile. Very thick lithospheres (∼100 km) would inhibit melt production to the extent that no outgassing would occur within 700 Myr.

Choice of initial volatile contents
Most of the bulk Earth's accreted volatile content may have been expelled during magma ocean degassing. The crystallised mantle's initial H 2 O and CO 2 contents, χ ini H 2 O and χ ini CO 2 in weight percent, are the amounts of volatiles remaining after degassing adopted from the magma ocean model of Elkins-Tanton (2008). The most water-rich scenario would correspond to just below an Earth ocean mass. These estimates are consistent with more recent modelling work by Barth et al. (0). Hier-Majumder and Hirschmann (2017) found comparable but slightly elevated values considering faster freezing of the magma ocean. Elkins-Tanton (2008) considers primordially-accreted volatile contents of 0.05-0.5 weight percent H 2 O and 0.01-0.1 weight percent CO 2 , although the true values are poorly constrained. Figure 1: The evolution of the extrusive volumetric melt production rate, assuming a melt density of 3000 kg m −3 . Note that these values represent only 10-40% of the total melt production. Cases are binned to 5%-increments in the extrusive volcanism fraction, and averaged over all other input parameters (Table 2). Lighter colours indicate higher extrusive volcanism fractions. Overlain in solid grey lines and corresponding to the secondary y-axis are the evolutions of the mean temperature over the whole mantle, for each individual model run.

Choice of mantle oxidation states
Geochemical clues mostly agree that since around 3.8 Ga or earlier, the Archean upper mantle was about as oxidised as present (Canil, 1997;Delano, 2001;Trail et al., 2011;Nicklas et al., 2018;Armstrong et al., 2019), although other analyses suggest a more gradual oxidation up to QFM (Aulbach and Stagno, 2016;Nicklas et al., 2019). Trail et al. (2011) investigated Hadean zircons with ages between 4.0-4.36 Ga and derived an oxidation state around QFM. If these zircons did indeed come from primitive mantle melts as suggested by Trail et al. (2011), this would suggest an oxidised mantle almost directly after the Moon-forming impact.
We allow mantle f O 2 to vary between the minimum Hadean value of three log-units below the equilibrium with the IW buffer (Wood et al., 2006), and a typical modern value of four log-units above the equilibrium. We assume a constant oxidation state, and later test the effect of log(f O 2 ) increasing linearly.

Choice of extrusive volcanism percentages
Only a fraction of magma will rise high enough to extrude and outgas at the surface. The remaining melt in the intrusive component will eventually replenish the mantle. Crisp (1984) estimate an extrusive volcanism percentage of 6-25% for intracontinental and island arc settings, and 16-33% at mid-ocean ridges and hotspots. Because this value is unknown for early Earth, we allow it to vary between 10% and 40% (Grott et al., 2011). We do not consider any diffuse degassing from intrusive volcanism, although fluid pressure would indeed transport intrusive magmatic gas to the surface through rock fissures.

Results
Before speciation, the total C (or H) outgassing flux is proportional to the product of the C (or H) melt concentration and the extrusive melt production rate. In principle, either of these factors has the potential to scale outgassing dramatically, although they do not necessarily show the same variance. We will first present the melt production rates and melt concentrations from the model. Then we will describe the resulting outgassing fluxes, their final speciation modulated by f O 2 .

Thermal history and melt production
The thermal histories underlying the outgassing model are shown in figure 1. Mean mantle temperatures gradually increase due to powerful primordial radiogenic heating and relatively sluggish stagnant lid convection. The temperature dependence of viscosity acts as a thermostat, such that by 700 Myr, average mantle temperatures converge to ∼2320-2400 K.
The first ∼50 Myr see spikes in the melt flux. Directly after magma ocean solidification, the mantle temperature is still very close to the solidus. In these conditions, slight mass movement from convection can trigger immediate re-melting. This depletes the residual mantle and raises the solidus temperature. Melting is revived once either radiogenic heating warms the upper mantle, or convection upwells undepleted material from the lower mantle to pressures of ≤12 GPa.
Note that higher initial temperatures than considered here would lead to unphysical conditions, with most of the mantle above the solidus temperature. Whilst a much wetter mantle might aid melting, we do not find a dramatic effect of χ ini H 2 O on the overall melt production. We do not include dehydration's stiffening effect on viscosity.

Melt contents of H 2 O and CO 2
Due to the different chemistry controlling the partitioning of C and H 2 O into the melt (section 2.3), χ melt H 2 O and χ melt CO 2 respond to different model parameters. These behaviours are reflected in figure 2. H 2 O partitioning into the melt depends directly on the local mantle H 2 O abundance, the maximum of which is set by χ ini H 2 O . The effect of mantle content on melt content is non-linear: in equation (7), χ melt H 2 O depends on the local melt fraction, whilst the presence of water facilitates melting by suppressing the solidus.
Carbon partitioning, in contrast, is strongly redox-dependent. This is true for the regime we model, where graphite is the stable phase, but would not apply to more oxidising conditions (Stagno, 2019). Note, that the local mantle source CO 2 concentration can still limit χ melt CO 2 . This effect is most obvious at very low values of χ ini CO 2 below 50 ppm, or for cases approaching maximum outgassing-in either situation the local carbon inventories can run out. The most carbon-rich simulations can be seen (figure 2) to reach peak χ melt CO 2 before 700 Myr for this reason. Meanwhile, the effect of f O 2 on χ melt CO 2 is drastic everywhere. With each step of one log-unit below the IW buffer, we see χ melt CO 2 drop by an order of magnitude. There is an early transient stage where the melt abundance drops steeply. This is associated with the early pulse of melting and hence depletion (figure 1).

Exploration of outgassing scenarios
A notable result of this work is that all scenarios, regardless of oxidation state, produce outgassing no higher than total estimates for modern Earth (e.g., Catling and Kasting, 2017). To see why this is, here we examine what controls outgassing rates in this model.

Outgassing as a function of mantle oxidation state
As expected (Holland, 1984;Kasting et al., 1993), there is a heavy dependence of the outgassing rate on f O 2 . Figure 3 shows the evolution of each species' mean outgassing flux in mol yr −1 , with cases binned by f O 2 and all other input parameters varying according to a random uniform distribution (Table 2). A more oxidised mantle is associated with more CO 2 and H 2 O outgassing, and a more reduced mantle is associated with more H 2 outgassing. The pattern for CO is more complicated. It peaks around IW: low f O 2 limits the amount of total carbonate that can dissolve into the melt, despite the gas-phase equilibrium favouring CO over CO 2 . Outgassing rates increase slightly with time (by 10-15%) during the 700 Myr modelled. The early pulse of outgassing is associated with the early pulse of melting (figure 1).
To illustrate the dependence of evolved outgassing on f O 2 , figure 4 shows every run's cumulative outgassed mass per species, as a function f O 2 . Both C-and H-bearing species are well-separated in terms of mass, except for cases close to the IW buffer. Increasing f O 2 from IW − 3 to IW + 4 is associated with an increase of over an order of magnitude in the mass of H 2 O, and an increase of at least five orders of magnitude for CO 2 , whilst the mass of H 2 decreases. The hydrogen speciation shows a gradient in redox. The carbon speciation also shows a gradient in redox, whilst the total amount of carbon depends on partitioning (itself a function of redox).
It is worth emphasising that figure 4 does not show the actual gas masses residing in the atmosphere. For instance, the atmospheric partial pressure of H 2 O is limited to the saturation vapour pressure of liquid water; on the Myr-timescales of this model, the cooling timescale of the outgassed plume to ambient atmospheric temperatures is irrelevant. At T surf = 333 K, the saturation vapour pressure is 0.2 bar. The cumulative mass of H 2 O (which would quickly condense) is at most about a tenth of modern Earth's ocean mass, which would support the idea that earlier magma ocean degassing supplied a significant fraction of the original ocean mass (Pahlevan et al., 2019).
The influence of f O 2 on the outgassed volatile masses is built into the model. Intrinsically, there are three reasons for the mass dependence in figure 4. The first two reasons comprise the redox-dependent speciation (setting total carbon) and volatile speciation (setting [H 2 O]/[H 2 ] and [CO 2 ]/[CO]) already described. The third reason is that the reduced gases investigated here have a lower molecular mass than their oxidised counterparts. Thus, for equal quantities outgassed, reduced atmospheres are necessarily less massive than oxidised atmospheres.       Figure 6 summarises the Spearman's rank (nonlinear) correlation coefficients between each species' outgassing flux, averaged over the final 10 Myr, and the eight input parameters plus the melt volumes and volatile concentrations. This section explains how the remaining correlations stem from the model.

Other factors influencing outgassing
Influence of initial temperatures. We find virtually no correlation between initial temperature conditions and outgassing fluxes after 700 Myr of convection. However, an initially hotter mantle closer to the solidus temperature forces a sharp pulse of melting in the first 100 Myr, expelling massive amounts of water and indeed raising the total mass of water outgassed. This is also not to say that instantaneous local temperatures are unimportant, as they affect the melt fraction as well as the equilibrium constants in (3), (10), and (11).
Influence of extrusive volcanism fractions. From (13), we expect an approximately linear relationship between the maximum outgassing flux and the fraction of melt allowed to contribute to outgassing. Although the correlation between this fraction and the computed outgassing fluxes is not as strong as f O 2 , accounting for extrusive volcanism in this model automatically downscales outgassing. For example, all else being equal, the difference between χ extr = 10% and χ extr = 100% would be an order of magnitude in outgassing.
Role of melt production rate versus melt volatile content. Figure 6 also shows the correlations of all outgassing fluxes with the instantaneous volatile melt concentrations and melt production rates. CO 2 outgassing fluxes correlate more strongly with CO 2 melt concentrations than with the total melt production rates, whilst the opposite appears to be true for H 2 O. Note that H 2 O melt concentrations are not independent of melt volumes because of partitioning of water into the melt-H 2 O concentrations are largest for the smallest melt fractions. To some degree, the relative strengths of the correlations reflect the experimental variances in melt production and volatile melt contents, given the ranges of the prior distributions on the input parameters. Nevertheless, in our simulations, outgassing will be ultimately limited by the volatile budget of the mantle. Thus rates of outgassing would not scale indefinitely with rates of volcanism. Figure 5 presents the experimental distributions of each species' cumulative outgassed partial pressures, shown for different redox bins. Again, the dominant effect of f O 2 is visible in how these distributions are separated along the x-axis. The median and 1σ values of these partial pressures are quantified in Table 3, along with the corresponding masses and fluxes.   and H-species outgassing is due to the mutual effect of f O 2 on both quantities, while the correlation of f O 2 with the total flux appears low because it does not affect the sum of H 2 and H 2 O.

Secular mantle oxidation
Several mechanisms for mantle oxidation have been proposed (e.g., Wood et al., 2006;Sharp et al., 2013;Wordsworth et al., 2018;Schaefer and Elkins-Tanton, 2018;Nicklas et al., 2019), including degassing itself (Kasting et al., 1993). These mechanisms would be associated with different durations. By tracking the concentration of abundant multivalent cations, one might self-consistently evolve the mantle redox state within our outgassing scenarios. We leave this to future work, and as a first step, simply consider a linear increase of log(f O 2 ) from IW − 2 to IW + 2 over 700 Myr, for a new set of 110 runs with otherwise random parameters as before (Table 2). Figure 7 demonstrates that secularly increasing log(f O 2 ) results in cumulative outgassing not obviously distinguishable from a constant log(f O 2 ) − IW ∈ [0, 2]. The spread in these cases would be largely explained by variations in the melting rate.
Low f O 2 does not affect melt concentrations of H-species, such that the same total amount of total H can be outgassed. Carbon, meanwhile, is strongly limited in melts during the early reduced stage, which might imply lower cumulative partial pressures. However, Figure 7 suggests that this effect is somewhat muted due to the scatter induced by other unknown parameters. Note, regardless, that small differences in cumulative outgassing could represent dramatic differences in the ultimate atmospheric composition; scenarios of evolving mantle redox likely deserve more rigorous treatment than attempted here.

Viscosity treatment
This work has considered a single set of upper mantle rheological parameters, corresponding to the canonical Arrhenius viscosity law for wet olivine from Karato and Wu (1993). Lower mantle rheology is taken from Tackley et al. (2013). Different viscosity treatments could potentially influence melting rates and therefore outgassing rates. However, several regulation mechanisms exist such that an immediate weakening of upper mantle viscosity, whilst increasing the convective velocity, may not lead to dramatically higher melting in the long term. Firstly, latent heat loss associated with any melting that does occur would cool the mantle (the latent heat of melting is indeed modelled here). These lowered temperatures would both limit further melting and stiffen the local viscosity in a negative feedback loop (Ogawa and Yanagisawa, 2011). Secondly, a sudden increase in convective velocity would lead to faster replenishment of depleted material in the upper mantle, which promotes melting; yet, thirdly, this effect would be offset simultaneously by more efficient cooling of the mantle, which suppresses melting.
The Karato and Wu (1993) law for wet diffusion creep does have a relatively weak temperature dependence and a strong pressure dependence, which could reduce some of the feedbacks described here. Nevertheless, as in Dorn et al. (2018), we expect that melting is ultimately controlled by the lower mantle viscosity because rapid flow in the upper mantle merely leads to rapid depletion. With that said, we would still emphasise that the results here are subject to our chosen viscosity treatment, and should be interpreted as such.

Submarine versus subaerial outgassing
We have assumed that all outgassing is subaerial and occurs at 1 bar. Several independent constraints place an upper limit on either the surface barometric pressure or pN 2 at 1.1 bar for 3.5-2.7 Ga (Marty et al., 2013;. A lack of plate motion might allow volcanic eruptions to construct topography above sea level relatively quickly (e.g., Olympus Mons, a Martian shield volcano comparable to Germany in area). In this case, volatiles would degas freely from surface lavas at atmospheric pressure. Because near 1 bar there is almost no effect on gas solubility or f O 2 in our chosen parameterisations, we do not consider solubility or any redox change during degassing itself. If most outgassing occurred on the seafloor, then the higher pressures would allow melts to retain more volatiles, as long as the gas and melt are in equilibrium (e.g., eruptions are not explosive). Global outgassing rates would then be lower than predicted here, and show a different ratio of outgassing species (Gaillard et al., 2011).

Note on CH 4 outgassing
We neglect the outgassing of CH 4 entirely because its vapour phase is not stable at the pressure, temperature, and redox ranges we consider (Zhang and Duan, 2009;Wetzel et al., 2013;Ramirez et al., 2014). Some CH 4 fluid could be produced from a reduced source at high pressures up to 11 GPa (Scott et al., 2004). Although CH 4 can also be stable at surface pressures, melt temperatures would be too high for it to degas directly.

Buoyancy limit of melting
As in Noack et al. (2014Noack et al. ( , 2017, we have not considered the possibility of melting at depths greater than 12 GPa. This assumes that basaltic melt is not buoyant above this pressure because it is denser than olivine (Ohtani et al., 1995). However, Mosenfelder et al. (2009) have suggested that basaltic melt can nevertheless remain less dense than ringwoodite and wadsleite up to ∼25 GPa. Therefore the actual density crossover pressure could be lower than modelled here (e.g., Beuchert and Schmeling, 2013). This could increase the volume of the mantle annulus in which melting can occur, possibly leading to enhanced outgassing if all buoyant melt at depth reaches the base of the lithosphere. However, it may be that no additional melting occurs at these depths, depending on how the local mantle temperatures compare to the (strongly pressure-dependent) melting temperature.

Stagnant lid outgassing effects
Some general consequences of imposing a stagnant lid regime are listed below. This compilation is by no means complete and would benefit from future work. Noack et al. (2014), Tosi et al. (2017), Foley and Smye (2018), and Gaillard et al. (2021) also discuss how tectonic regimes can affect outgassing rates. The arguments here could be tentatively supported by the observation that Venus currently appears less volcanically-active than Earth (Smrekar et al., 2010). show the 1σ deviation. For context, we include estimates of modern Earth's CO 2 outgassing, spanned by blue rectangles for the mid-ocean ridge (MOR) system, and by beige rectangles for hotspots. MOR estimates are taken from the data in Le Voyer et al. (2019) and , where the ranges of outgassing fluxes and melt production rates correspond to their quoted uncertainty on the global total, and the ranges of mantle and melt contents correspond to the 2σ-width of their log-normal sample distribution. For hotspots, the CO 2 flux lower limit is the sum from major hotspots from , and the generous upper limit is taken from Marty and Tolstikhin (1998); ocean island basalt (OIB) mantle source and melt concentrations are from ; magma supply rates are the sum over 19 hotspots from Mjelde et al. (2010). Estimates of the mantle source f O 2 are from O'Neill et al. (2018) and Amundsen and Neumann (1992) for MORB and OIB respectively.
i. Interior temperatures: For a constant mass, stagnant lid planets have hotter mantles because they lack the efficient cooling by subducting plates. This would enhance melting and volcanism. ii. Melt volumes: However, stagnant lid planets may be associated with less continuous melting for the same interior temperature because, three-fold, the melt zone is separated from the surface by a thick lithosphere, convection is weaker due to the lower temperature contrast, and the depleted mantle is not refertilised by subducting plates. iii. Volatile inventories of the magma source: Higher mantle volatile contents would lead to higher melt contents insofar as they reflect each other. The mantle source supplying stagnant lid volcanism should be less volatile-rich than the source of arc volcanism, which contains subducted crustal material relatively heavy in carbonates, organic carbon, and water (Wallace, 2005). The mid-ocean ridge source is drier, yet still generally appears carbon-rich compared to the initial concentrations we have adopted for the early degassed, differentiated mantle . Hotspots, meanwhile, sample the deep mantle, which seems to contain both lingering nebular and recycled components (e.g., Miller et al., 2019). Analyses of ocean island basalts suggest even higher mantle source concentrations than mid-ocean ridges-likewise associated with high-CO 2 melts (e.g., Shorttle et al., 2015;Miller et al., 2019), as are ridge segments in proximity to hotspots (Le Voyer et al., 2019). iv. Upper mantle hydration: If subducted oceanic crust is wet, then the presence of water would depress the solidus (Katz et al., 2003). This could facilitate melting near subduction zones in a plate tectonics regime.

Modern Earth's volcanic outgassing
A key result from this work is that our stagnant lid outgassing model produces average outgassing rates no higher than those estimated for modern Earth under plate tectonics. Whereas ∼68% of our cases more oxidised than IW + 3 have a CO 2 flux of between 0.6 and 4.1 Tmol yr −1 , estimates for Earth's global CO 2 outgassing now rest at around 8.5 ± 2 Tmol yr −1 (Catling and Kasting, 2017). Metamorphic degassing  Dasgupta and Hirschmann (2010) for CO 2 hotspot volcanism, multiplied by 7-8, the typical ratio of H 2 O to CO 2 in hotspot volcanic gas from Holland (1984); i.e., following Catling and Kasting (2017). Mantle source concentrations refer to Hawaii estimates from Wallace (1998); melt concentrations are also Hawaii estimates from Wallace and Anderson (1998). Mid-ocean ridge volcanism is omitted because it is known to be very dry, compared to the range shown here. Other data are explained in the caption to figure 8. and arc volcanism, irrelevant to a stagnant lid regime, both account for a quarter of this global estimate. Although our model does not necessarily apply to hotspot or mid-ocean ridge outgassing in a plate tectonics regime either, estimates of the modern rates of these processes might provide a useful point of context for illustrating the mechanisms that keep our outgassing results relatively low.
Namely, we see four main factors potentially influencing outgassing: mantle volatile content, melt production rates, melt volatile contents, and mantle oxidation state. To this end, figures 8 and 9 summarise our modelled outgassing rates as a function of these four factors, overplotted by the relevant present-day estimates for hotspot and ridge volcanism.
It should be noted upfront that our model may not be able to correctly capture modern ridge and hotspot outgassing rates, especially for CO 2 . As explained in section 4.1.5 above, these processes sample recycled material generally having both (i) higher oxidation states, at the QFM buffer or above (Amundsen and Neumann, 1992;Moussallam et al., 2016;O'Neill et al., 2018;Moussallam et al., 2019), to which our carbon partitioning model does not apply, and (ii) higher source volatile concentrations than our quasi-primordial bulk mantle values . The combined effect of (i) and (ii) is that CO 2 can potentially reach much higher melt concentrations on average. In oxidised conditions where carbonate is stable, abundant carbon would be lost very efficiently during partial melting.
Further to this point, modern outgassing rates are quite difficult to measure. Most quoted outgassing rates come from multiplying estimates of the melting rate with estimates of the concentration of volatiles in the melt. The volatile contents of melts are largely uncertain because melts are partially degassed; estimates rely on the behaviour of geochemical proxies (see Michael and Graham, 2015). An in-depth review of Earth's ridge and hotspot outgassing rates is well beyond the scope of this work.
Mantle oxidation state. Whilst shown here for completeness, mantle f O 2 does not play an important role in controlling modern CO 2 or H 2 O outgassing because the modern upper mantle is oxidised essentially everywhere, such that H 2 and CO would never be favoured, and carbon partitioning in the melt is not redox-dependent.
Melt production. Estimates of globally-integrated magma supply rates at mid-ocean ridges tend to show better agreement-most not far from 20 km 3 yr −1 (Crisp, 1984;Mjelde et al., 2010;Le Voyer et al., 2019)although these estimates refer to the volume of all melt contributing to crust production, which may not necessarily be the same as the volume of melt contributing to outgassing. This study has found stagnant lid extrusive melting rates up to tenfold higher. Therefore melting rates alone would not explain our model's low outgassing fluxes. Hotspots are more geographically sparse, and their global magma supply rate is lower than ridges.
Mantle concentrations. The ridge mantle source shows orders-of-magnitude variation in CO 2 concentrations (as derived from melt concentrations). This variation, extending from a minimum of 10 ppm to a maximum of 1980 ppm , encompasses our entire range of χ ini CO 2 . The hotspot mantle source tends to be even more volatile-rich, up to well over 1000 ppm in CO 2 , and several hundred ppm for H 2 O (Wallace, 1998;Miller et al., 2019).
Melt concentrations. Our oxidised cases overlap fairly well with estimates of the typical volatile concentrations in modern ridge and hotspot magmas, as well as with the corresponding outgassing fluxes. Simulations that reach higher χ melt CO 2 and χ melt H 2 O , with respect to modern Earth, accordingly tend to produce higher CO 2 and H 2 O outgassing rates. However, for many of our simulations, the mantle is almost completely depleted in volatiles by 700 Myr. Thus even if the convective velocity is increased by a reduction in viscosity, for example, we would not expect outgassing rates to increase any further.
In summary, the most reliable way to attain higher outgassing rates in our stagnant lid model would be to raise the volatile content of the melt. However, a volatile-rich melt is not necessarily attainable by tweaking model parameters because the melt partitioning of CO 2 and H 2 O will be ultimately limited by their stocks in the upper mantle source. To illustrate what a stagnant lid regime's lower interior volatile supplies mean for possible outgassing rates: a mantle with χ ini CO 2 = 180 ppm and no return fluxes could outgas a maximum of ∼16 bar CO 2 at 30% depletion and χ extr = 40%. Meanwhile, outgassing 8.5 Tmol yr −1 CO 2 (Catling and Kasting, 2017) over 1 Gyr is equivalent to a cumulative ∼75 bar.

Archean outgassing proxies
Few observational constraints exist on Archean outgassing rates. Xe isotope anomalies in Archean quartz suggest that mantle degassing was about tenfold greater at 3.3 Ga than at present (Avice et al., 2017;Marty et al., 2019). This proxy would apply to C-O-H outgassing rates at 3.3 Ga if they were derived directly from magma production rates-indeed, these high melting rates suggested by Xe isotopes are matched by our stagnant lid model. However, the carbon and water contents of the melt are at least as important as the volume of the melt in determining C-O-H outgassing rates.

From outgassing to atmospheric composition
The atmospheric composition and oxidation state do not mimic the volcanic gas composition or oxidation state: photochemical reactions and H escape can oxidise reduced gases, and extraterrestrial impactors could provide significant reducing power . In addition, it is possible that an earlier atmosphere degassed from the magma ocean was in some part still present at the end of the Hadean (Hamano et al., 2013;Nikolaou et al., 2019;Stüeken et al., 2020). We do not model the atmospheric composition here because its complexities deserve more astute attention. However, to build our discussion of this work's potential implications for the Archean Earth system, we will briefly contextualise our results alongside previous work. We focus on processes affecting the partial pressure of CO 2 due to its importance in long-term climate stability.
The equilibrium between outgassing sources and silicate weathering sinks (seafloor and subaerial) controls pCO 2 . Weathering is thought to regulate the surface temperature via a negative feedback loop: in the classic Walker et al. (1981) framework, there would be a single equilibrium combination of pCO 2 and surface temperature for fixed values of the outgassing flux and other relevant parameters (e.g., the weatherable surface area). Lower outgassing-all else being equal-implies a cooler climate, a weaker need to weather, and a lower CO 2 weathering flux (e.g., Kadoya and Tajika, 2014;Krissansen-Totton et al., 2018). For example, Krissansen-Totton et al. (2018) find pCO 2 between 2-500 mbar for a contemporaneous outgassing rate of 3-9 Tmol yr −1 (95% confidence intervals). According to that model, pCO 2 < 2 mbar might therefore be anticipated for less than 3 Tmol yr −1 of outgassing; such rates are not ruled out by this study. Note that this pCO 2 estimate assumes that weathering is not limited by mineral supply, and following equation (29) in Foley and Smye (2018), the supply limit to silicate weathering would only be reached at CO 2 fluxes well over 100 Tmol yr −1 in stagnant lid scenarios (i.e., it is not reached). For comparison, proxy analyses by Hessler et al. (2004) give a lower limit of pCO 2 ∼ 2.5 mbar at at 3.2 Ga.

Greenhouse warming under the Faint Young Sun
Despite the fainter luminosity of the young sun, geochemical evidence exists for a temperate climate and stable oceans on Earth as early as 4.4 Ga (e.g., Wilde et al., 2001;Valley et al., 2002Valley et al., , 2014. Attempts to explain this paradox often invoke stronger-than-modern atmospheric partial pressures of greenhouse gases (see review in Charnay et al., 2020). Volcanic outgassing, being the main source of these greenhouse gases, is a primary control on steady-state pCO 2 and surface temperature (e.g., Walker et al., 1981;Sleep and Zahnle, 2001;Kadoya and Tajika, 2014;Höning et al., 2019).
Indeed, the CO 2 outgassing rates predicted here are at best several times lower than the values often assumed in Archean Earth climate models (e.g. Sleep and Zahnle, 2001;Wordsworth and Pierrehumbert, 2013;Charnay et al., 2017;Kanzaki and Murakami, 2018;Krissansen-Totton et al., 2018). These higher rates might relate to some combination of a different tectonic mode at the time, a higher ratio of extrusive to intrusive melt, and a more-oxidised mantle.
3D general circulation models can be used to estimate the minimum partial pressures of CO 2 that would provide the necessary greenhouse warming, for a given solar luminosity. For example, Wolf and Toon (2014) find that maintaining 15°C requires ∼200 mbar of CO 2 at 3.8 Ga, if other climate parameters (e.g., day length) are not optimal, and tens of mbar if they are. Note that this smaller value is similar to the minimum pCO 2 from the low-outgassing scenarios studied by Krissansen-Totton et al. (2018). Dedicated atmospheric modelling should estimate the feasibility of greenhouse warming fed by weaker outgassing in the early Archean.
At the least, if a planetary mantle f O 2 were below the IW buffer, then it is apparent that virtually no CO 2 could be outgassed in principle, even with melting rates much higher than estimated here. Indeed, silicate Earth may have already reached f O 2 > IW + 1 by the end of the Hadean (Pahlevan et al., 2019). Depending on the longevity of the primordial magma ocean atmosphere, requiring a CO 2 outgassing rate of at least several Tmol yr −1 to sustain an early greenhouse may lend a constraint on the timing of mantle oxidation.

Conclusions
This work has coupled a 2D numerical model of stagnant lid convection with melting, volatile partitioning into the melt, and chemical speciation of these volatiles. The model returns redox-dependent volcanic outgassing fluxes of CO 2 , CO, H 2 O, and H 2 . We have applied this model to Hadean-Archean Earth. We find global CO 2 and H 2 O outgassing fluxes on the order of 1 Tmol yr −1 and 10 Tmol yr −1 respectively, depending on the mantle oxidation state. These fluxes are kept low by the assumption of a stagnant lid regime, wherein outgassing may never be much stronger than predicted here if the upper mantle cannot be efficiently replenished in volatiles. Our model may not apply to a scenario where most outgassing occurred on the seafloor, or where the upper mantle were too oxidised for graphite to be the stable form of carbon.
Coupled convection-outgassing models could inform further studies of the early Earth's atmospheric evolution. In particular, unknown CO 2 outgassing rates throughout the late Hadean and early Archean account for a significant portion of the uncertainty on the contemporaneous greenhouse warming capacity. Previous estimates of Earth's climate state during the early period of lower solar luminosity tend to employ outgassing rates higher than the present day, an assumption which might not be substantiated given the unknown tectonic state. If we believe that the presence of surface liquid water in the early Archean demands fairly high CO 2 partial pressures (Charnay et al., 2020), then coupled outgassing models might test whether an early initiation of plate tectonics-combined with rapid mantle oxidation-may be key to building a temperate young planet.