Dynamics of the Great Oxidation Event from a 3D photochemical-climate model

From the Archean toward the Proterozoic, the Earth's atmosphere underwent a major shift from anoxic to oxic conditions, around 2.4 to 2.1 Gyr, known as the Great Oxidation Event (GOE). This rapid transition may be related to an atmospheric instability caused by the formation of the ozone layer. Previous works were all based on 1D photochemical models. Here, we revisit the GOE with a 3D photochemical-climate model to investigate the possible impact of the atmospheric circulation and the coupling between the climate and the dynamics of the oxidation. We show that the diurnal, seasonal and transport variations do not bring significant changes compared to 1D models. Nevertheless, we highlight a temperature dependence for atmospheric photochemical losses. A cooling during the late Archean could then have favored the triggering of the oxygenation. In addition, we show that the Huronian glaciations, which took place during the GOE, could have introduced a fluctuation in the evolution of the oxygen level. Finally, we show that the oxygen overshoot which is expected to have occurred just after the GOE, was likely accompanied by a methane overshoot. Such high methane concentrations could have had climatic consequences and could have played a role in the dynamics of the Huronian glaciations.

The best constraints on the GOE come from sulfur isotopic ratios in precambrian rocks (Farquhar et al., 2007;Lyons et al., 2014). In the Archean anoxic atmosphere, the sulfur photochemistry was responsible for mass-independent fractionation of sulphur isotopes in sedimentary rocks (Farquhar et al., 2000).
The loss of mass-independent fractionation in sedimentary rocks less than 2.5 Gyr ago is explained by a modification of sulfur photochemistry main pathways due to the rise of the amount of oxygen, starting from 10 −5 present atmospheric level (PAL) (Pavlov and Kasting, 2002). This is directly linked to the amount of UV flux received for the sulfur photochemistry, which could be reduced by the appearance of ozone at higher oxygen level .
The GOE represents a major event in the history of the Earth. It profoundly impacted the atmospheric and oceanic chemistry, the climate, the mineralogy and the evolution of life. O 2 was a poison for a lot of anoxygenic form of life supposed already developed. Consequently, the GOE likely induced a retreat for anoxygenic form of life, including heterotrophic methanogens (i.e. organisms producing methane). Methane may have been more abundant in the anoxic Archean atmosphere than today (1.8 ppm), with levels as high as 10 4 pmm according to Catling and Zahnle (2020). Such high methane concentrations would have produced a significant greenhouse effect. A variation of 10 times the abundance of methane changes the mean surface temperature of about 4 K according to Charnay et al. (2020). Furthermore, Sauterey et al. (2020) show that the diminution of CH4 combined with the regulation of CO2 by the carbonate-silicate cycle favors the triggering of an ice age. The decrease of the biological methane productivity and the methane photochemical lifetime could have reduced its abundance and thus its warming contribution, potentially triggering the Huronian glaciations that took place between 2.4 and 2.1 Gyr (Kasting et al., 1983;Haqq-Misra et al., 2008) The GOE is a key event to understand the co-evolution of life and environment on Earth but also on exoplanets. However, major questions remain concerning the triggering and the dynamics of the GOE.
Before the appearance of oxygenic photosynthesis, the redox state of the atmosphere was ruled by the balance between reductant fluxes from volcanism and metamorphism and the hydrogen atmospheric escape (Catling et al., 2001). The appearance of oxygenic photosynthesis, which was much more efficient than previous mechanisms of photosynthesis, relying on ubiquitous H 2 O, CO 2 and light, profoundly changed the biogeochemical cycles. The oxygen is produced by oxygenic photosynthesis (summarized by formula (1)), which can be reversed by aerobic respiration.
The burial of organic carbon (a very small fraction of the net primary productivity) allows the accumulation of oxygen until that an equilibrium is reached between the burial of organic carbon, the reductant fluxes, the oxidative weathering (i.e. the oxidation of buried organic carbon re-exposed to the surface by plate tectonics) and the hydrogen escape. Assuming a methane-rich atmosphere, atmospheric oxygen is also strongly coupled to methane through the reaction of methane oxidation: On the early Earth and once oxygenic photosynthesis appeared, methane would have been mostly produced by the fermentation of organic matter followed by acetogenic methanogenesis: For aerobic conditions, it could have been consumed by oxygenic methanotrophy, similar to reaction (2). Goldblatt et al. (2006) and Claire et al. (2006) developed simplified models of the biogeochemical cycles of O 2 and CH 4 . They proposed scenarios for the evolution of the amount of oxygen and methane and the dynamics of the GOE.
Several hypotheses have been proposed to explain the mechanisms and timing of the GOE. The hydrogen escape is one of them, proposed by Catling et al. (2001). Considering the irreversible oxidation of methane by the reaction (4) : and the reverse reaction (2), we get a chain of reaction that causes a net gain of oxygen by transforming 2H 2 O into O 2 (5) : Emergence of continents and subaeriel volcanism is another hypothesis developed in Gaillard et al. (2011) that led to a change of the chemical composition and the oxidation state of sulfur volcanic gases, precipitating the atmospheric oxygenation.
But whatever precipitated the GOE, the rise of oxygen seems to be linked to an atmospheric instability caused by the formation of the ozone layer and its impact on the photochemical methane oxidation (Goldblatt et al., 2006). Slowly increasing O 2 , by oxygenic photosynthesis, would have accumulate enough starting the ozone layer formation. The ozone layer provide a photochemical shield which limit the oxygen photochemical destruction leading to the methane oxidation. Therefor the oxygen could have accumulate more easily, producing more ozone, shielding more efficiently the oxygen destruction and then rising an instability of growing oxygen until others processes would have limited the oxygen abundance, such as rock oxidation.
In this paper we focus on the atmospheric photochemical losses by methane oxidation associated to the GOE. Previous 1D studies of Goldblatt et al. (2006); Claire et al. (2006); Zahnle et al. (2006), developed dynamical models of the GOE based on 1D photochemical models. In this study, we use for the first time a 3D Global Climate Model to compute the chemical lifetime of the different species and to explore 3D effects. Since the oxygen build-up is linked to the formation of the ozone, we could expect effects from the latitudinal/longitudinal ozone distribution or from the variations of UV irradiation by the seasonal and the diurnal cycles. We also investigate the potential links between the GOE and the Huronian glaciations. In particular, the consequences of a cold climate (.e.g. a snowball Earth event) on the photochemistry have never been studied.
Following, we describe the atmospheric model used for this study in Section 2. In Section 3, we analyze the photochemistry of the late Archean/Neoproterozoic with the 3D model, highlighting the chemical lifetimes and the impact of the global mean surface temperature. Based on these results, we describe the dynamical evolution of O 2 and CH 4 during the GOE in Section 4, highlighting consequences of the Huronian glaciations. We finish with a summary and perspectives in Sect. 5.

Model
A 3D photochemical global climate model is used to characterize photochemical oxygen losses in the atmosphere, dominated by methane oxidation in the model. The model, the LMD-generic, is a generic Global Climate Model (GCM) initially developed at the Laboratoire de Météorologie Dynamique (LMD) for the study of a wide range of atmospheres. It allows to model easily different atmospheres, which makes it widely used, for instance to study early climates in the solar system D. Wordsworth et al., 2012;Charnay et al., 2013;Turbet et al., 2020b) or extrasolar planets (Selsis et al., 2011;Leconte et al., 2013;Bolmont et al., 2016;Fauchez et al., 2019;Turbet et al., 2020a). From the photochemical module of the martian version of the model (Lefèvre et al., 2004), we have developed a generic version. It allows to easily adapt the chemical network and it introduces the calculation of the photolysis rates and their heating rates within the model.
The chemical network is derived from the REPROBUS model of the present Earth stratosphere (Lefèvre et al., 1998) but adapted to the assumed composition. Halogen, heterogeneous and nitrogen chemistry is not taken into account due to the weak constraints available. Furthermore, halogen and heterogeneous chemistry have a negligible effect on the oxygen chemistry studied. In contrast, the chemical network includes a detailed methane chemistry. It allows to take into account the different pathways of the methane oxidation balance. The methane network is built according to Arney et al. (2016) and Pavlov et al. (2001). The whole chemical network is detailed in appendix A.
We have developed a new photochemical module for the LMD-generic code. Although previous versions of the code already included photochemical processes, they were hard-coded to specific atmospheres. The module is now flexible and no longer uses pre-computed photolysis rates, which are now computed using the actual absorbers abundances in the model. The module also accounts for the heating rates by photolysis although the abundance of O 3 in the present study is too low to yield significant heating. Including the heating by photodissociations will nevertheless be essential for other potential applications of the generic model, including Earth-like oxygen-rich atmospheres.
The GCM is adapted to the supposed conditions of the Archean Earth. The rotation period is set to 17.5 hrs (adapted for an orbital period of 500 days according to Zahnle and Walker (1987); Bartlett and Stevenson (2016) -1 day equal to 17.5 hrs ).
The spectrum of the star is calculated for 2.7 Ga and 1.0 AU from Claire et al. (2006) (see Figure 1). We define an atmosphere with 98% of N 2 and 1% of CO 2 for 1 bar at the surface. CH 4 has been added to the radiative contribution with different levels using the HITRAN 2012 database. The topography of the Archean Earth presumes a central continent. We define then a ocean planet with an equatorial supercontinent as in Charnay et al. (2013) The 1D version of the model uses a surface ocean, a surface mean albedo of 0.28, a mean solar zenith angle of 60 • and an eddy diffusion coefficient vertical profile from Zahnle et al. (2006).
Photochemical losses from the atmosphere by methane oxidation in the GCM are balanced with a production flux at the surface. This flux is established by fixing the abundance of the species considered in the first layer of the GCM (formula (8)). It gathers all the surface contribution, such as organic carbon burial or input of reductants, which can be in a dynamical equilibrium with atmospheric chemistry. When the stationary regime is reached, we can quantify the total photochemical loss/production of a species as its integrated surface flux required to maintain constant its surface abundance. Several simulations are performed on a grid of oxygen and methane abundance at the surface. The calculated fluxes are used to determine the photochemical loss flux as a function of the variable oxygen and methane abundances at the surface. Gregory et al. (2021) suggested to be cautious on the surface boundary conditions concerning a fixed mixing ratio or a fixed flux. While Gregory et al. (2021) does not describe the oxygen instability area with the fixed flux boundary condition, the fixed mixing ratio boundary condition allow us to explore the full equilibrium states. Fixing the flux is indeed more physical (or in this case biological) than fixing mixing ratios as flux ratio can vary in time and space during the evolution toward the steady state. But in the present study we wanted on purpose to determine the fluxes that correspond to a given steady state composition.
The resulting fluxes are finally consistent with a 1/2 ratio driven by the oxygen and methane chemistry. Beyond this, additional simulation done by fixing the flux have shown some differences, which can not fully be explained by the analysis of Gregory et al. (2021) considering the 3D geometry of the surface but the homogeneous boundary condition for each approach. This need to be properly analyze in a study focusing only on this aspect, and could be done following this paper.
The GCM ensures the convergence of carbonaceous species by also fixing the CO 2 abundance at the surface to 1%. In addition, the GCM takes into account atmospheric escape according to the Catling et al. (2001) model. The H 2 abundance in the last layer of the GCM is updated according to the formula 6 thanks to the escape flux calculated by the formula 7.

Methane oxidation fluxes
The atmospheric oxygen loss is dominated by the oxidation of methane, formula (2). The oxidation of methane is catalyzed by OH radicals (Gebauer et al., 2017). These radicals are produced by the photochemistry of water vapor: The amount of water vapor in the troposphere is controlled by temperature in a 1D model with an infinite water reservoir on the surface but in a 3D dynamic model with dry continental surfaces, it also depends on the horizontal transport, evaporation and precipitations.
Photochemical processes depend on insolation and therefore on diurnal and seasonal variations. The formation of the ozone layer is a turning point for the photochemical balance. The ozone layer produces an UV shield, which limits the photochemical processes leading to methane oxidation and destroying oxygen. Oxygen can accumulate more efficiently and form more ozone.
This positive feedback creates an oxygen instability which accounts for the sudden oxygenation of the atmosphere and may therefore be sensitive to the spatial distribution of ozone, and thus to the global 3D transport.
In this section, we compare the results between the 1D and 3D model on photochemical oxygen and methane losses at steady state. We compute the variation of vertical chemical pathways to methane oxidation as a function of the surface O 2 fixed abundance. We also analyze the spatial distribution of ozone. Finally, we examine how surface fluxes required to sustain a steady state depend on surface temperature.

From 1D to 3D models
We ran the 1D model until steady state for a range of O 2 from 10 −7 to 10 −3 volume mixing ratio and CH 4 from 10 −6 to 10 −3 volume mixing ratio . Figure 2 shows the total atmospheric O 2 loss (F O2 ) as a function of these two parameters. These results are consistent with the previous study of Zahnle et al. (2006).  Differences can come either from averaging the UV irradiation geometry in 1D or from horizontal and vertical transport. A priori, the photochemical losses (F O2 and F CH4 ) are not significantly modified ( Figure 3) and the vertical transport seems responsible for these differences. The 3D vertical transport seems to transport species more efficiently than the 1D transport model which uses an Eddy coefficient from Zahnle et al. (2006) to mimic the 3D transport. In particular, this results in smaller vertical gradients with the 3D model ( Figure 4).
Despite of the aforementioned small departures, we find that the 1D model reproduces the results of the more comprehensive 3D model.

Vertical distribution of O 2 losses
At steady state (averaging seasonal variations) O 2 photochemical losses compensate for the surface outgassing of O 2 . The ratio 1/2 between F O2 and F CH4 (Figure 3) shows that the atmospheric losses of oxygen and methane are dominated by the methane oxidation reaction 2, which uses two molecules of O 2 for each molecule of methane (and forms one molecule of CO 2 and 2 molecules of H 2 O). While oxygen molecules are mainly involved with CO and CO 2 cycle (Gebauer et al., 2017), CH 4 is dominated by the methane oxidation. O 2 , CO and CO 2 cycle is not a global loss of oxygen because there is no sources/losses of CO at the surface and consequently O 2 has a null balance, contrary to the methane oxidation. So, to analyze how losses of molecular oxygen are distributed in altitude, it is clearer to look at the methane loss, dominated by the oxidation of methane.
The Figure 5 represents the rate of methane destruction as a function of altitude and for different O 2 levels. As for their integrated value, the vertical profiles of F O2 and F CH4 are similar when computed with 1D and 3D models are similar. We distinguish the contribution to the loss of three main altitude domains : the troposphere, stratosphere and above. The losses are dominated, whatever the O 2 abundance, by the tropospheric contribution. In the appendix B, we identify empirically the main reaction pathways leading to a net methane oxidation 2. Figure 5 shows a migration of losses from the troposphere to the stratosphere when oxygen abundance increases. Catalysis by OH remains at the heart of the oxidation mechanism although a different reaction pathway is identified (Appendix B): as the stratospheric ozone becomes more abundant, the production of O( 1 D) through its photolysis increases: which then initiates the production of OH through the reaction 10 instead of the O( 1 D) coming from O 2 photolysis. The stratospheric contribution is less efficient than the tropospheric one, that is why F O2 decreases for the high abundances of   Figure 6 shows how the ozone column density during a year varies with latitudes from low to high pO 2 . At low pO 2 we observe a maximum of ozone column density close to the poles during winter, while it switch to summer and at lower latitudes (∼20 • ) at high pO 2 . For comparison, we computed the ozone column density for the current Earth with and without continent. Our simulation of the current Earth with continent reproduces quite well the present-day ozone distribution (Tian et al., 2010), with a maximum close to the north pole in March and at 50 • south in October. The pole-ward transport of ozone in our simulation seems weaker than in reality, certainly due to the absence of the effect of gravity waves, which play a major role in the Brewer-Dobson circulation. By removing continent, the ozone maximum occurs later for both hemispheres, mostly due to the thermal inertia of the ocean. For such a case, the maximum of ozone occurs at the same season as the early Earth with high pCO 2 but at higher latitudes (∼45 • ). This is due to the faster rotation of the early Earth, limiting the latitudinal extension of the stratospheric circulation.

Ozone
The 3D effects have no significant consequence on the photochemical losses of O 2 . The decrease in F O2 with increasing levels of O 2 above an O 2 vmr of 10 −5 (Figure 3) is due to the formation of the ozone layer, due to UV shielding of the methane oxidation (Goldblatt et al., 2006;Zahnle et al., 2006). Figure 6 shows latitudinal variation relatively less than a factor 2, which have a tiny impact on the ozone lifetime and the photochemistry, and figure 7 shows that ozone is relatively homogeneous over the whole planet making 1D modeling relevant. Nevertheless, there are some variations of the O 3 column with latitude and related variations of the biologically harmful UV flux reaching the surface (Figure 8). This non-homogeneous opacity of the ozone UV shield may be important for the evolution and distribution of organisms living at the surface.

Surface temperature effect
The resulting tropospheric temperature profile , calculated in both 1D and 3D models, follows a moist/dry adiabat that controls the vertical profile of water vapour, which, in turn, affects the greenhouse warming. Since this water vapor is the source of OH, catalyzer of methane oxidation, this interplay suggests a possible link between surface temperature and photochemical losses of CH 4 and O 2 , which we investigate here.
The previous results ( Figure 3) were obtained with a surface temperature close to 280 K, whatever the oxygen abundance at the surface and an abundance of methane of 10 −4 . The surface temperature is the same because all the parameters are the same (rotation period, obliquity, solar input, surface composition (water/continent), continental albedo, etc...) and main greenhouse gases (CO 2 , CH 4 ) are not changed depending on the oxygen level.
As a first test to assess the effect of surface temperature, we use the 1D model with a surface surface temperature forced to remain at 220 K (by setting to zero the heating/cooling rate of the surface). This is the approximate surface temperature that would be reached by the 1D model with a frozen start, although the actual value would depend on the level of greenhouse CH 4 . This way we can evaluate the impact of a snowball event on the photochemistry. Figure 9 shows that photochemical losses decrease with the temperature. How strong this decrease is depends on the O 2 level, the drop being the largest around the maximum of F O2 , so for a vmr of O 2 around 10 −5 . At this O 2 level, the photochemical losses occur entirely in the  troposphere. Figure 10 compares the vertical profiles of methane photochemical losses for a surface temperature of 280 K and 220 K and for different oxygen abundances at the surface. We see that the influence of surface temperature on the losses is located in the troposphere. Thus, the larger the tropospheric contribution, the larger the decrease in losses. The stratospheric thermal profile and the tropopause (the cold trap that controls the transport of water vapor to the stratosphere) show little response to a decrease of surface temperature from 280 K to 220 K (see Figure 13). As a consequence only the tropospheric chemistry is affected.
The main implication of this result is that, under conditions of a global glaciation, the oxygen instability is triggered for an oxygen abundance/flux about an order of magnitude lower. For a given O 2 flux, glacial conditions should favor the switch from an oxygen-poor to an oxygen-rich atmosphere. During glaciations, however, environments able to provide both light and   The temperature trend of oxygen losses is determined. These losses are calculated using the 1D model by setting different surface temperatures for a methane abundance of 10 −4 and an oxygen abundance of 10 −5 at the surface. Figure 12 shows F CH4 (equivalent to 1 2 F O2 ) as a function of surface temperature. We observe an order of magnitude increase in losses per 50 K. The discontinuity at 273 K is artificially produced by the change of albedo at the surface between ice-free oceans (albedo = Figure 11. Species profile for 1D model surface temperature 280K (solid) and 220K (dashed). 0.07) and fully ice-covered oceans (albedo = 0.65). The increase in albedo below 273 K reflects more UV flux and increases photochemical losses in the troposphere. The 3D model smooths this effect by gradually freezing the ocean at the surface.
We use the slow convergence of the 3D model towards a frozen state. The photochemical equilibrium is established on a time scale of a few years, whereas the progressive freezing of the surface takes place over several decades. A quasi-stationary state of species abundance in the atmosphere and consequently of atmospheric losses is rapidly established. During freezing on a longer time scale, the evolution of temperature and oxygen losses in their quasi-stationary state is recorded to establish the link between oxygen losses and temperature, see Figure 12. The 3D model converges to a frozen state where the sea ice extends from the poles to 20-25 • N/S. The coverage is about 60-65%, as observed in Charnay et al. (2013), and the surface temperature converges to 230 K. A cold start with a completely frozen surface is then performed to evaluate the impact of the 3D model is weaker than the 1D model. At the root of this difference is the difference in transport between the two models. The 3D model mixes the tropospheric temperature more efficiently, warming the troposphere and reducing the impact of a cooling on atmospheric losses, see Figure 13.  Temperature appears to have a significant impact on atmospheric loss Between a temperate state and an ice age, oxygen losses can vary by a factor of 2 to 4. Oxygen abundance and methane abundance can accumulate more in the atmosphere during an ice age. The consequences of ice ages on the temporal evolution of oxygen and methane during the GOE are studied in the following by considering this new trend.
After the GOE, a carbon-13 isotopic variation of nearly 15 ‰ is observed (Lyons et al., 2014). This event is called the Lomagundi event (Bachan and Kump, 2015). Although the dynamics of the oxygenation process remain uncertain, this event suggests an over-oxygenation of the atmosphere (Catling and Zahnle, 2020). Other evidence supports this phenomenon, such as the evolution of the δ 34 S fraction of carbonate-associated sulfates (Planavsky et al., 2012;Schröder et al., 2008) or fluctuations in the degree of uranium enrichment in organic-rich shales (Partin et al., 2013). An increase in overall oxygen productivity followed by a decrease would seem to account for this over-oxygenation (Harada et al., 2015;Holland and Bekker, 2012;Hodgskiss et al., 2019). Nevertheless, we have seen previously that a global glaciation phenomenon can significantly decrease the atmospheric oxygen losses. It is therefore not impossible that Huronian glaciations could also have had an impact on atmospheric over-oxygenation during the global oxygenation process. In order to study the dynamics of over-oxygenation, a time evolution equation model is established based on the equations of Goldblatt et al. (2006), Claire et al. (2006) and adjusted to take into account variations in primary oxygen productivity and atmospheric losses. We compare the previous Goldblatt et al. (2006) parametrisation for atmospheric loss with the GCM interpolated function on the time evolution without an over-oxygenation using the time evolution equation model established. We then apply this model to reproduce the over-oxygenation and the fluctuations brought by glaciations.

Equations model
The temporal evolution of oxygen and methane abundance during the GOE is modeled by the Goldblatt et al. (2006) equations.
They relate atmospheric losses, described in section 3, to surface contributions associated with biogeochemical exchanges. Goldblatt et al. (2006) propose for this biosphere model a parameterization of atmospheric losses according to Ψ O2 [CH 4 ] 0.7 where Ψ O2 = 10 a1ψ 4 +a2ψ 3 +a3ψ 2 +a4ψ+a5 , ψ = log([O 2 ]), a 1 = 0.0030, a 2 = −0.1655, a 3 = 3.2305, a 4 = −25.8343 and a 5 = 71.5398. Goldblatt et al. (2006) also estimate the value of the different parameters (see Table 1). A set of values is associated with a steady state of oxygen and methane abundances. Depending on the value of the flux of reductant r, one is on a state of oxygen-poor or oxygen-rich equilibrium. In order to study the temporal evolution between these two equilibrium states, we use the results of Claire et al. (2006) to establish a temporal evolution of the flux of reductant r (Figure 14). We introduce this temporal evolution in the equations of Goldblatt et al. (2006) to reproduce the dynamics of oxygenation. Finally, we introduce a coefficient α N (≥ 1) to model a photosynthetic over-productivity responsible for the over-oxygenation, as well as a coefficient α Ψ (≤ 1) to model the decrease of the atmospheric losses during a glaciation. The evolution of the abundance of oxygen, methane and buried carbon is then described by equations: where the different terms and values are detailed in Table 1 The size and shape of the over-oxygenation is uncertain, as are the possible variations in oxygen sources and losses during this process. We then arbitrarily define a log-polynomial fit for the time evolution of the parameters α i = e (a(t−t0) 2p +b) (α i = α N or α Ψ ). When there is no deviation from the initial model the values of α i are equal to 1. The variation of the α i parameters is triggered in a way that is consistent with the predictions of previous studies on the evolution of primary productivity and ice ages.

Dynamics with constant temperature and constant primary productivity
The Goldblatt et al. (2006) parameterization for atmospheric oxygen loss as a function of oxygen and methane abundance is compared (Figure 15) to the results obtained with the 1D GCM, which are shown to be similar to the 3D model. The Figure 14. Temporal evolution of the reductant contributions (see Table 1) with αi = 1 and using for oxygen atmospheric loss (FO 2 ) the Goldblatt et al. (2006) parametrization.
parameterization established by Goldblatt et al. (2006) does not seem to correctly capture the methane dependence. This variation is not independent of oxygen abundance and cannot be described by a constant power x of methane abundance At high oxygen abundances (> 10 −4 ) the variation seems to increase with increasing methane, and conversely at lower oxygen abundances (< 10 −4 ). An asymptote appears to emerge at low oxygen abundances for the highest methane abundances (> 10 −5 ). This behaviour is already seen in Zahnle et al. (2006). and GCM 1D results (dash). Figure 16 shows the equilibrium states as a function of the total reductant parameter r with the Goldblatt et al. (2006) parameterization and an interpolation of the 1D GCM results for the atmospheric losses. These curves show the stable and unstable equilibrium states of the atmosphere. They justify the rapid switch from an oxygen-poor to an oxygen-rich state. The flux of reductant that triggers the instability varies from about 1. × 10 10 to 3. × 10 10 mol O 2 equiv. yr −1 with the interpolation of the 1D GCM results. The asymptotic behaviour of the atmospheric losses calculated with the 1D GCM at low oxygen and high methane (high r) levels, can be seen on the equilibrium states with a approximatevly constant value at those levels. This behaviour remains uncertain. A time evolution of oxygen and methane abundance with these two atmospheric loss models ( Figure 17) shows that with the interpolation of 1D GCM results oxygenation is faster. The equilibrium positions Figure 16 show that indeed the amount of oxygen is more sensitive to the flux of reductant. Oxygenation therefore occurs for a smaller variation of the reductant flux. As the methane abundance is directly related to the reductant flux ([CH 4 ] = r s ), we also observe a smaller variation of the methane abundance during oxygenation.

Long overshoot with variable primary productivity
We model a long variation over 400 million years of the primary productivity, parameter α N , with a first phase of overproduction before a return to the initial production. We represent Figure 18 the reference evolution for constant α N equal to 1, as well as two different intensities of over-productivity with a maximum at 2 and 10 times the initial productivity. The reference model is in good agreement with the results of Goldblatt et al. (2006) and Claire et al. (2006). We observe a variable over-oxygenation depending on the intensity of the primary productivity variation, but also an over-abundance of methane.
The primary productivity corresponds to a photosynthetic production of oxygen but also of organic matter, transformed then into methane by the methanogens. Consequently, the production of methane is increased as well as that of oxygen. Such an over-abundance of methane is not highlighted by previous studies, nor by the geological record. It is difficult to constrain the amount of methane at that time. This scenario is therefore possible, although one could also imagine that oxygen enrichment of the atmosphere limited the conversion of organic matter into methane. A conversion that is carried out by methanogens, developing more favorably in a reducing environment. Similar explanation has been proposed for the debated higher level of methane during the Boring Billion periods, 1.8-0.8 Ga (Pavlov et al., 2003). This temporal evolution is identical if we model an inverse variation of atmospheric losses (α Ψ = 1 α N ), corresponding to a glaciation event of 400 million years. This scenario is less likely since there are several shorter glaciation events, called Huronian glaciations (Young et al., 2001). In addition, a change in primary productivity provides a link to the positive anomaly in carbon-13 isotope fractionation (Lyons et al., 2014). By using the model of the evolution of the isotopic ratio of Goldblatt et al. (2006), we can evaluate the impact of the variation of primary productivity or atmospheric losses on this ratio: where f is the fraction of buried volcanic carbon, δ carb is the carbon-13 isotope ratio δ 13 C for carbonates, δ i is the carbon-13 isotope ratio δ 13 C for volcanic carbon, and δ org is the carbon-13 isotope ratio δ 13 C for organic carbon. The fraction of buried volcanic carbon can be described by the biosphere model based on the fraction C of buried carbon where C ∼ β w (N + r). We establish the approximation that f α (N + r). Goldblatt et al. (2006) defines the initial state during the Archean with δ carb = 0 ‰, δ i = -6 ‰ and δ org = -30 ‰ giving f = 0.2. With this initial state, we determine the relation of proportionality between f and (N + r), then we evaluate the evolution of the isotopic fractionation of the carbonates δ carb during the temporal evolution of the oxygen and methane abundance thanks to the formula: Figure 19 represents the evolution of δ carb as a function of time for the time evolution models established with a primary over-productivity reaching a factor of 2 and the inverse evolution of atmospheric losses. A primary over-productivity is consistent with the anomaly of about 15 ‰ measured in contrast to the decrease in atmospheric losses. Figure 19. Temporal evolution of δ carb . Models with αN reaching a factor 2 and αΨ reaching a factor 0.5.

Short overshoot with variable temperature
The Huronian glaciations represent several glacial events that took place during the oxygenation of the atmosphere at the beginning of the Proterozoic. The model of the α Ψ parameter variation over 400 million years is therefore not consistent with the episodic nature of these glaciations. To reflect the impact of these glaciations on atmospheric losses, an episodic variation over shorter times during the oxygenation period can be established. Figure 20 presents the temporal evolution of oxygen and methane following this adjustment. We observe the punctual over-oxygenations linked to the variations of α Ψ with a global trend that follows the temporal evolution of the initial model where α i = 1. We note, during the over-oxygenation, an increase in the methane abundance. This increase in greenhouse gases can trigger the thawing of the surface. During the thawing, the surface temperature increases causing an increase in atmospheric losses by oxidation of methane. This is why the amount of methane decreases again. This negative feedback, coupled with the hysteresis phenomenon between the frozen and warm state of the surface, could be one of the key factors to explain the cyclic character of Huronian glaciations.  Figure 21 shows the complete model for over-oxygenation with a 400 million year variation in primary productivity and a maximum factor of 2, as well as a fluctuation provided by shorter time scale variations in atmospheric losses due to Huronian glaciations with a maximum factor of 0.5. This coupling is not intended to exactly reproduce what happened during the GOE but to highlight both effects on the evolution of oxygen and methane during this event.

Summary
The atmospheric equilibrium states during the Archean and the Proterozoic has been established for the first time using a 3D photochemical-climate model. Despite some 3D transport discrepancies the atmospheric-surface equilibrium fluxes of methane and oxygen are not significantly different from a 1D model as it has been done in Zahnle et al. (2006) or Gebauer et al. (2017).
Following, the (photo)chemical equilibrium pathways have been determined depending on the altitude. It highlights the evolution of the tropospheric and stratospheric contribution depending of the oxygen levels. What remains however constant is the important link with the OH molecules which catalyze the methane oxidation. Because of that, we found a crucial dependency of the surface fluxes equilibrium with the surface temperature. This temperature dependence is sensitive to the 3D transport and appears weaker in 3D than in 1D. However, a global glaciation could reduce the oxygen atmospheric losses by a factor 2 to 10. Taking this new contribution into account in a time evolution model, we show that glaciations bring fluctuations in oxygen and methane abundance with an overshoot during glaciations. The increase in methane following glaciation produces an additional greenhouse effect that could eventually lead to deglaciation. This warming increases again the atmospheric losses of methane and it is possible to establish a cycle of glaciation-deglaciation. Few evidences of methane (Lowe and Tice, 2004;Cadeau et al., 2020) have been yet discovered but it could explain how the planet has been warmed up enough to terminate the glaciation thanks to the methane greenhouse effect. The feedback between the glaciation and the methane evolution coupled with the glaciation hysteresis process might also explain the multiple glaciations known as the Huronian glaciations. More generally, the link between temperature and photochemical processes shows that a decrease in temperature favors the oxygenation of the atmosphere. Without falling into a global glaciation, phenomena such as the decrease of the atmospheric CO 2 or the emergence of continents induce a decrease of temperature favorable to the oxygenation of the atmosphere.

Perspectives
Beyond the temperature, the UV-Visible stellar radiation received by the planet controls the photochemical processes and the quantity of ozone, essential for the oxygenation phenomenon of the atmosphere. A red dwarf such as TRAPPIST-1a presents a spectral distribution in favor of the UV radiation with respect to the visible radiation which favors the accumulation of ozone.
The 1D study of Gebauer et al. (2018) shows that around such red dwarfs the oxygenation of the atmosphere is triggered for lower atmospheric oxygen levels and surface oxygen flux. These irradiation conditions then favor the oxygenation of the atmosphere. A compact system such as TRAPPIST-1 has presumably synchronous planets (Vinson and Hansen, 2017), which brings an interest to use 3D models. The permanent dichotomy between day and night side brings an important variation of temperature  and photochemical processes. The red dwarfs represent the majority of the stars and would thus shelter the majority of the planets. It is then necessary to understand how an oxygenated atmosphere could evolve around these planets. For this, 3D models will be necessary to capture the effect of a synchronization. . If k1a,0, k 1b,0 , k1a,∞ and Fc specified, the formula is k1a,0 Ogawa and Ogawa (1975) O2 and  and Gibson et al. (1983) and Minschwaner et al. (1992) and Yoshino et al. (1988) and Fally et al. (2000) CO2 Chan et al. (1993) CO2 and Stark et al. (2007) and Yoshino et al. (1996) and Parkinson et al. (2003) and  O3 Sander et al. (2006) O3 et al. (2005) and Chung et al. (2001) and Thompson et al. (1963)  Competing interests. The authors declare that they have no conflict of interest.