© ESO 2018

1 Introduction

Planetary atmospheres evolve due to interactions with the planet’s surface and losses into space. At the surface, gas can be removed from the atmosphere by several processes, such as subduction (Marty & Dauphas 2003), and added to the atmosphere by other processes, such as outgassing during magma ocean solidification (Noack et al. 2014). At the top of the atmosphere, gases are lost to space, which over time can lead to significant atmospheric erosion (Lammer et al. 2014; Luger et al. 2015). Atmospheric loss into space takes place by a large number of different mechanisms (e.g. Lammer et al. 2008). One factor that is common to almost all of these processes is the fact that the loss rates depend strongly on the thermal and chemical structure of the upper atmosphere. Atmospheres that are hotter and more expanded have higher loss rates by essentially all mechanisms (e.g. Lichtenegger et al. 2010).

Much recent work has studied hydrodynamic losses of atmospheres. Many of these studies concentrate mostly on atmospheres composed primarily of H and He (e.g. Lammer et al. 2014; Shaikhislamov et al. 2014; Luger et al. 2015; Khodachenko et al. 2015; Owen & Mohanty 2016). Such atmospheres can experience very high hydrodynamic losses, largely due to the small average molecular masses of the gas (Erkaev et al. 2013). Atmospheres composed of water vapour, such as the possible early atmosphere of Venus, are also likely to undergo hydrodynamic escape as the dissociation of H₂O creates large amounts of atomic H (Lichtenegger et al. 2016). Generally more interesting for planetary habitability are atmospheres dominated by heavier molecules, such as CO₂, N₂, and O₂. The physical processes in these atmospheres are very complex (e.g. Kulikov et al. 2007; Tian et al. 2008a), and detailed models are needed to understand their structures. Such atmospheres are less likely to undergo hydrodynamic losses due to their higher molecular masses, and other atmospheric loss processes must be taken into account, such as polar ion outflows (Glocer et al. 2012; Airapetian et al. 2017) and pick-up of exospheric gas by the stellar wind (Kislyakova et al. 2014) and coronal mass ejections (Khodachenko et al. 2007; Lammer et al. 2007). In all cases, the specific atmospheric composition is critically important for the detailed physics of the upper atmosphere (Kulikov et al. 2007).

The most important input into the upper atmospheres of planets is the irradiation by the central star, especially in X-ray and ultraviolet (together XUV)¹ wavelengths, though IR photons can also be important. The absorption causes dissociation and ionization, and significant heating. The energy gained by this heating is mostly lost by cooling due to IR emission from several molecules, most notably CO₂. In the upper thermosphere of the Earth, the local heating is much stronger than the local cooling, and the excess energy is transported into the lower thermosphere by thermal conduction. The chemical structure of the upper atmosphere is determined by the composition of the lower atmosphere, chemical/photochemical reactions, and diffusion. Sophisticated models that take into account all of these processes have been applied for solar system planets for decades (e.g. Fox & Bougher 1991; Roble 1995; Ridley et al. 2006), but only a few studies have applied such models to planetary atmospheresunder very different conditions to those of the current solar system terrestrial planets (Tian et al. 2008a; Tian 2009).

The need for sophisticated first principles upper atmosphere models is clear when considering the range of atmospheric conditions that exist. In addition to different atmospheric compositions, the distribution of planets spans the entire range of possible masses and orbital distances from their host stars (López-Morales et al. 2016). Furthermore, different planets are exposed to very different conditions from the central star. Observations of young solar analogues have shown that the Sun was much more active in X-rays and UV than it currently is (Güdel et al. 1997; Ribas et al. 2005). Recently, Tu et al. (2015) showed that the early evolution of the Sun’s activity depended sensitively on its early rotation rate; this is important since we do not know how rapidly the Sun was rotating, and different evolutionary tracks for XUV can lead to different atmospheric evolution scenarios (Johnstone et al. 2015b). Stellar activity evolution depends also on the star’s mass, with lower mass stars remaining highly active for longer amounts of time (West et al. 2008). In addition, the exact shape of a star’s XUV spectrum depends on its spectral type and activity (Telleschi et al. 2005; Johnstone & Güdel 2015; Fontenla et al. 2016).

The aim of this paper is to develop and validate a first principles physical model for the upper atmospheres of planets and to apply it to the Earth to understand how the atmosphere reacts to changes in the CO₂ abundances and the solar XUV spectrum. This physical model will be used as an important component in future studies on the evolution of terrestrial atmospheres. In Sect. 2, we present the complete physical model. In Sect. 3, we validate the model by calculating the atmospheric structures of Earth and Venus. In Sect. 4, we study the effects of enhanced CO₂ abundances and the effects of the solar XUV evolution between 3 Gyr in the past and 2.5 Gyr in the future on the structure of the Earth’s upper atmosphere. In Sect. 5, we summarise and discuss our results. To solve the physical model presented in this paper, we have developed The Kompot Code, which we describe in the Appendices². In the appendices, we describe in detail the numerical methods used to solve the equations described in Sect. 2.

2 Model

2.1 Model overview

The purpose of our model is to calculate the atmospheric properties as a function of altitude for arbitrary planetary atmospheres. The input parameters are the planetary mass and radius, the atmospheric properties at the base of the simulation, and the stellar radiation spectrum at the top of the atmosphere. Our computational domain is 1D and points radially outwards from the planet’s centre, extending between the lower boundary at an arbitrary altitude in the middle atmosphere to the upper boundary at the exobase. In the description of the state of the atmosphere, we make a few basic assumptions. Firstly, we assume that the gas has one bulk advection speed shared by the entire gas, though different chemical species have different diffusion speeds. Secondly, we assume that the neutrals, ions, and electrons have their own temperatures that evolve separately. Thirdly, we assume quasineutrality, meaning that the electron density is equal to the total ion density everywhere. The two stellar inputs are the XUV (i.e. X-ray and ultraviolet) field between 10 and 4000 Å, and the infrared field between 1 and 20 μm.

In this model, we break the gas down into components in two separate ways: in Eq. (1) the gas is broken down by different chemical species (e.g. N₂, O₂, CO₂, etc.), and inEqs. (3)–(5), the gas is broken down into neutrals, ions, and electrons. In the rest of the paper, we define the “components” of the gas as the neutral, ion, and electron gases, and are referred to using the subscripts n, i, and e. Unless otherwise stated, when we discuss the electrons, we are referring to the thermal electron gas, and not the non-thermal electrons produced in photoionization reactions.

The main physical processes taken into account in this model are

• atmospheric expansion/contraction in response to changes in the gas temperature and composition,

• the transfer of X-ray, ultraviolet, and infrared radiation through the atmosphere, including the production of non-thermal electrons by photoionization reactions,

• atmospheric chemistry, including photochemistry and reactions driven by impacts with non-thermal electrons,

• molecular and eddy diffusion,

• neutral heating by stellar XUV and IR radiation,

• electron heating by impacts with non-thermal electrons,

• infrared cooling, particularly by CO₂ molecules,

• heat conduction for each gas component,

• and energy exchange between the components.

The equations that describe the changes of the atmosphere due to these processes are $$\begin{array}{lll}\frac{\partial n_j}{\partial t}+\frac{1}{r^2}\frac{\partial \left[r^2(n_{j}v+{\Phi }_{\text{d},j})\right]}{\partial r}\hfill & =\hfill & S_j,\hfill \\ \frac{\partial (\rho v)}{\partial t}+\frac{1}{r^2}\frac{\partial \left[r^2\left(\rho v^2+p\right)\right]}{\partial r}\hfill & =\hfill & -\rho g+\frac{2p}{r},\hfill \\ \frac{\partial e_{\text{n}}}{\partial t}+\frac{1}{r^2}\frac{\partial \left[r^{2}v\left(e_{\text{n}}+p_{\text{n}}\right)\right]}{\partial r}\hfill & =\hfill & -{\rho }_{\text{n}}vg\left(Q_{\text{h},\text{n}}-Q_{\text{c},\text{n}}-Q_{\text{in}}-Q_{\text{en}}\right)\hfill \\ \hfill & \hfill & +\frac{1}{r^2}\frac{\partial }{\partial r}+[r^2{\kappa }_{\text{mol}}\frac{\partial T_{\text{n}}}{\partial r}\hfill \\ \hfill & \hfill & +r^2{\kappa }_{\text{eddy}}\left(\frac{\partial T_{\text{n}}}{\partial r}+\frac{g}{c_{\text{P}}}\right)],\hfill \\ \frac{\partial e_{\text{i}}}{\partial t}+\frac{1}{r^2}\frac{\partial \left[r^{2}v\left(e_{\text{i}}+p_{\text{i}}\right)\right]}{\partial r}\hfill & =\hfill & -{\rho }_{\text{i}}vg+\left(Q_{\text{h},\text{i}}-Q_{\text{c},\text{i}}-Q_{\text{ei}}+Q_{\text{in}}\right)\hfill \\ \hfill & \hfill & +\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2{\kappa }_{\text{i}}\frac{\partial T_{\text{i}}}{\partial r}\right],\hfill \\ \frac{\partial e_{\text{e}}}{\partial t}+\frac{1}{r^2}\frac{\partial \left[r^{2}v\left(e_{\text{e}}+p_{\text{e}}\right)\right]}{\partial r}\hfill & =\hfill & -{\rho }_{\text{e}}vg+\left(Q_{\text{h},\text{e}}-Q_{\text{c},\text{e}}+Q_{\text{ei}}+Q_{\text{en}}\right)\hfill \\ \hfill & \hfill & +\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2{\kappa }_{\text{e}}\frac{\partial T_{\text{e}}}{\partial r}\right],\hfill \end{array}$$

where r is the radius, n_j is the number density of the jth species, ρ is the total mass density, v is the bulk advection speed, ρv is the momentum density, ρ_k, e_k, p_k, and T_k are the mass density, energy density, thermal pressure, and temperature of the kth component of the gas, Φ_d,j and S_j are the diffusive particle flux and chemical source term of the jth species, g is the gravitational acceleration, Q_h,k and Q_c,k are the heating and cooling functions for the kth component, Q_ei, Q_in, and Q_en are the electron–ion, ion–neutral, and electron–neutral heat exchange functions, κ_mol and κ_eddy are the molecular and eddy thermal conductivities, κ_i and κ_e are the ion and electron thermal conductivities, and c_P is the specific heat at constant pressure. Since chemistry and diffusion do not change the total mass density of the gas, Eq. (1) implies the standard mass continuity equation. It is also important at times to calculate γ (i.e. the ratio of specific heats) for the neutral and ion gases, which are mixtures of species with different γ values; for this we assume γ_j = 5∕3 for atomic species and γ_j = 7∕5 for molecular species³.

The exobaseis assumed to be where the mean-free-path of particles becomes larger than the pressure scale height. The mean free path is calculated from l_mfp = 1∕(σN), where σ is the total collision cross-section, and N is the total number density. In reality, different species have different σ; however, the values tend to be similar and our calculated exobase location is not sensitive to small changes in σ. We therefore assume σ = 2 × 10⁻¹⁵ cm⁻² always.

2.2 Hydrodynamics

Including gravity, the purely hydrodynamic parts of Eqs. (1)–(5) are $$\begin{array}{lll}\hfill & \hfill & \frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial \left(r^2\rho v\right)}{\partial r}=0,\hfill \\ \hfill & \hfill & \frac{\partial (\rho v)}{\partial t}+\frac{1}{r^2}\frac{\partial \left[r^2\left(\rho v^2+p\right)\right]}{\partial r}=-\rho g+\frac{2p}{r},\hfill \\ \hfill & \hfill & \frac{\partial e_{\text{n}}}{\partial t}+\frac{1}{r^2}\frac{\partial \left[r^{2}u\left(e_{\text{n}}+p_{\text{n}}\right)\right]}{\partial r}=-{\rho }_{\text{n}}ug,\hfill \end{array}$$

where the ion and electron energy equations are identical to the neutral energy equations with the n subscript replaced with the i and e subscripts. In Appendix B, we give an explicit method for solving these equations. Explicitly solving the full set of hydrodynamic equations is undesirable when the atmosphere is static, or close to static. We note that no atmosphere is ever fully hydrostatic since there is always some escape at the top of the atmosphere, meaning that there will always be a net upward flow of material.

In this paper, we use a method for solving the hydrodynamic equations given by Tian et al. (2008a). This method is not appropriate when the atmosphere is transonic, which is often the case for strongly irradiated planets. The basic simplifying assumption of the method is that the mass and momentum density structures are in a steady state, such that ∂ρ∕∂t = 0 and ∂(ρv)∕∂t = 0. Eqs. (6) and (7) can then be written $$\begin{array}{lll}\hfill & \hfill & 2r\rho v+r^{2}v\frac{\text{d}\rho }{\text{d}r}+r^2\rho \frac{\text{d}v}{\text{d}r}=0,\hfill \\ \hfill & \hfill & \frac{2\rho v^2}{r}+2v\rho \frac{\text{d}v}{\text{d}r}+v^2\frac{\text{d}\rho }{\text{d}r}+\frac{\text{d}p}{\text{d}r}=-\rho g.\hfill \end{array}$$

Assuming an ideal gas, the pressure is given by $p={\overline{m}}^{-1}\rho k_{\text{B}}T$ and the radial derivative of p is $$\frac{\text{d}p}{\text{d}r}=-\frac{\rho }{\overline{m}}{v}_0^2\frac{\text{d}\overline{m}}{\text{d}r}+v_0^2\frac{\text{d}\rho }{\text{d}r}+\rho \frac{v_0^2}{T}\frac{\text{d}T}{\text{d}r},$$(11)

where $v_0^2=k_{\text{B}}T/\overline{m}$ and $\overline{m}$ and T are the average molecular mass and temperature of the entire gas. Putting these three equations together gives $$\begin{array}{lll}\hfill & \hfill & \frac{1}{v}\left(1-\frac{v^2}{v_0^2}\right)\frac{\text{d}v}{\text{d}r}=\frac{1}{T}\frac{\text{d}T}{\text{d}r}+\frac{g}{v_0^2}-\frac{1}{\overline{m}}\frac{\text{d}\overline{m}}{\text{d}r}-\frac{2}{r},\hfill \\ \hfill & \hfill & \frac{1}{\rho }\frac{\text{d}\rho }{\text{d}r}=-\frac{1}{T}\frac{\text{d}T}{\text{d}r}-\frac{g}{v_0^2}+\frac{1}{\overline{m}}\frac{\text{d}\overline{m}}{\text{d}r}-\frac{v}{v_0^2}\frac{\text{d}v}{\text{d}r}.\hfill \end{array}$$

When v = 0, Eq. (13) gives the density structure of a hydrostatic atmosphere. As in Tian et al. (2008a), we solve these equations consecutively. Firstly, we update the energies using Eq. (8). With the new temperature structure and the already known structure of $\overline{m}$, we then recalculate the structure of v by integrating from the exobase downwards through the grid using Eq. (12), assuming the outflow speed at the exobase is known. Finally, since the density at the lower boundary of the simulation is a fixed value and is therefore known, we calculate the structure of ρ by integrating upwards to the exobase using Eq. (13). For both the solution of the energy equations and the integration of ρ and v, we use the implicit Crank–Nicolson scheme, as described in Appendix D.

When the atmosphere is supersonic at the upper boundary, the material has escape velocity and simply flows away from the planet; in these cases, the appropriate boundary conditions are zero-gradient outflow conditions. Specifically, the values for each quantity in the final grid cell are made equal to the values in the second to last grid cell. When the atmosphere is subsonic at the upper boundary, we assume an outflow speed that is consistent with the Jeans escape rate. We first calculate the Jeans mass escape rate, Ṁ_Jeans, using the expressions given in Luger et al. (2015), and then calculate the upper boundary velocity from $v_{\text{exo}}={\dot{M}}_{\text{Jeans}}/(4\pi r_{\text{exo}}^2{\rho }_{\text{exo}})$.

When the bulk flow of the atmosphere is not negligible, the effects of advection on the species densities must be taken into account. We do this using the advection scheme described in Appendix B by converting the calculated cell boundary mass fluxes into individual species particle fluxes. When using the semi-static hydrodynamic approach, we use the advection scheme to calculate the mass fluxes only (i.e. Eqs. (B.11)–(B.18)). Since the total mass density structure is being calculated at every timestep assuming that it has already come to a steady state, the changes in the density are a result of advection that is not explicitly calculated in the model. To take this into account, after updating the structure of ρ using Eq. (13), we scale the species number densities by a species-independent factor at each location to ensure that ρ = ∑_jm_jn_j, where the sum is over all species.

2.3 Stellar radiation and non-thermal electrons

2.3.1 X-ray and ultraviolet radiation

The most important external input into the upper atmosphere is the star’s X-ray and ultraviolet (=XUV) radiation. Its importance stems from the fact that atmospheric gases absorb radiation at XUV wavelengths very effectively. This means that the XUV radiation is absorbed high in the atmosphere where the gas densities are low and relatively small energy inputs can lead to large temperature changes. The XUV spectrum also drives the most important chemical processes in the upper atmosphere, and is therefore essential for calculating the chemical structure of the atmosphere.

We irradiate the atmosphere with a stellar XUV spectrum between 10 and 4000 Å. The XUV spectrum is divided into 1000 energy bins and represented by the irradiance, I_ν, which is theenergy flux per unit frequency. This input spectrum is assumed to be unattenuated at the exobase. The radiation transfer through the atmosphere is then calculated based on the density structures of each absorbing species.

A weakness of 1D atmosphere models is that in reality the planet is being irradiated from one side only, which makes fully simulating the atmosphere at minimum a 2D problem. In 1D models, approximate simplifying assumptions must be made. We assume that the computational domain is pointing in an arbitrary direction relative to the position of the star. The angle between this direction and the direction that points directly at the star is the zenith angle, θ. We calculate the XUV spectrum at each point in the atmosphere by doing the radiation transfer from the exobase to each point separately. This geometry is demonstrating in Fig. 1, where the dashed black line shows the path that the radiation takes through the atmosphere. We assume that the state of the atmosphere at any given altitude is uniform over all latitudes and longitudes. This means that when doing the radiation transfer, we get the densities of each species at each given point by taking the values at the point in our simulation that has the same altitude. In all simulations in this paper, except our Venus simulation in Sect. 3.2, we assume a zenith angle of 66°. We find in Sect. 3.1 for the case of the Earth that this gives a decent representation of the atmosphere averaged over all longitudes and latitudes.

To calculate the XUV spectrum at a given grid cell, we integrate along the dashed black line in Fig. 1 for each energy bin using spatial stepswith length Δs given by H∕5, where H = N∕(dN∕dr) is the density scale height. The change in the irradiance over a path Δs is given by $$I_{\nu }(s+\Delta s)=I_{\nu }(s){\text{e}}^{-\Delta {\tau }_{\nu }},$$(14)

where Δτ_ν is the optical depth along the path length Δs and is given by $$\Delta {\tau }_{\nu }=\Delta s{\displaystyle \sum _j{\sigma }_{\nu ,j}}[R_j],$$(15)

where the sum is over all photoreactions in our chemical network, σ_ν,j is the cross-section of the jth photoreaction at frequencyν, and [R_j] is the number density of the reactant in the jth photoreaction. By summing over individual photoreactions instead of using the total absorption cross-sections for each species, we ensure that the radiation transfer and the photochemistry are fully consistent. The individual photoreaction cross-sections are discussed in Sect. 2.4.1. In Fig. 2, we show the XUV spectrum at several altitudes in our model for the current Earth.

Fig. 1 Cartoon illustrating the geometry of the radiation transfer in our model. The bottom right of the cartoon is the centre of the planet and the incoming stellar radiation is travelling horizontally from left to right. The inner black region is the planet and the shaded region is the upper atmosphere. The simulation domain of our 1D atmosphere model, illustrated with the blue region, is a region that is pointing radially outwards from the planet centre. The angle θ is the zenithangle. To calculate the stellar XUV and IR spectra at any given point, we perform radiation transfer along the dashed black line.

Fig. 2 Stellar XUV irradiance spectrum at different altitudes in the atmosphere of the Earth, calculated using our current Earth model presented in Sect. 3.1. The black and purple lines show the spectrum at the top and bottom of our model respectively.

2.3.2 Infrared radiation

Another input into the atmosphere that can be important is the stellar infrared radiation. Although this has a negligible effect on the Earth’s upper atmosphere, it is a signficant source of heating for Mars and Venus (Bougher & Dickinson 1988; Fox & Bougher 1991). This difference is due to the different abundances of CO₂, which is a strong absorber of IR radiation. We calculate the transfer of the stellar IR spectrum between 1 and 20 μm through the atmosphere, and its effect on atmospheric heating. For the input stellar IR spectrum, we assume a simple black-body spectrum with a temperature of 5777 K. We make the same geometrical assumptions for the IR radiation transfer as we make for the XUV, as demonstrated in Fig. 1, and perform the integration from the exobase to each grid cell using the same method described for XUV. For IR transfer, the sum in Eq. (15) for the optical depth is over all considered absorbing species, where σ_ν,j and [R_j] are the cross-sections and number densities of the jth absorbing species. We consider only absorption of IR radiation by CO₂ and H₂O molecules. In future studies, the influences of other molecules will be included when necessary.

We calculate the absorption spectra of CO₂ and H₂O using the software package kspectrum (Eymet et al. 2016), which is an open-source code for calculating the high resolution absorption spectra of common atmospheric gases using the HITRAN 2008 and HITEMP 2010 molecular spectroscopic databases (Rothman et al. 2009; Rothman et al. 2010). Although the absorption spectrum is temperature and pressure dependent, we calculate the cross-sections at 200 K and 10⁻² mbar only and use these values everywhere in the atmosphere. In order to resolve all features in the CO₂ absorption spectrum, a large number (~10⁶) of spectral bins are needed. The wavelength-dependent absorption cross-sections for CO₂ and H₂ O are shown in Fig. 3. However, including so many energy bins is computationally too expensive for our model; instead, we only consider energy bins that have cross-sections above 10⁻²² cm² for at least one of the considered molecules. Tests have shown that we get identical results using this threshold, while limiting the number of energy bins to something reasonable (~ 10⁴).

The heating of the atmosphere by the absorption of IR radiation is discussed in Sect. 2.5.1. We calculate the heating assuming that all energy removed from the radiation field by absorbtion is immediately added to the thermal energy reservoir of the neutral gas. What actually happens is that the absorption of photons excites the molecules and the heating of the gas only takes place when they are then collisionally deexcited. Some of this energy will not in fact end up in heat, but will be reradiated back into space. To take this into account, we add an additional excitation term into the equations for 15 μm CO₂ cooling in Sect. 2.5.2. We write the excitation rate due to stellar IR photons as $$S_{\text{IR}}=\frac{{\displaystyle {\int }_{\nu }{\sigma }_{\nu ,{\text{co}}_2}}[{\text{CO}}_2]I_{\nu }\text{d}\nu }{{(h\nu )}_{15\mu \text{m}}},$$(16)

where the integral is over all considered frequencies, ${\sigma }_{\nu ,{\text{co}}_2}$ and [CO₂] are the absorption cross-section and number density of CO₂ at frequency ν, and ${(h\nu )}_{15\mu \text{m}}$ is the energy of a 15 μm photon. The cgs units for S_IR are excitations s⁻¹ cm⁻³. The numerator in Eq. (16) gives the volumetric heating rate due to the absorbtion of IR photons by CO₂. The assumption here is that all energy absorbed from the IR field by CO₂ eventually contributes to the excitation ofthe 15 μm bending mode in CO₂ molecules. The main absorption bands are at 15, 4.3, 2.7, and 2.0 μm. The latter two are combination bands, and photons absorbed in these bands cause multiple excitations in the 15 and 4.3 μm band transitions. Our assumption in Eq. (16) is reasonable if the majority of energy in the 4.3 μm vibrational state is transferred to the 15 μm vibrational state by vibrational–vibrational exchanges, as argued by Taylor & Bitterman (1969; see the discussion in Sect. 2 of Dickinson 1976).

Fig. 3 Infrared absorption spectra of CO2 (upper panel) and H2O (lower panel) as calculated by kspectrum (Eymet et al. 2016). The gas is assumed to have a temperature of 200 K and a pressure of 10−2 mbar. The horizontal dashed line shows the cross-section of 10−22 cm−3.

2.3.3 Non-thermal electrons

Many of thephotoionization reactions that take place in the upper atmosphere are caused by photons that contain significantly more energythan is needed simply to cause the ionization. This additional energy is given to the produced photoelectrons in the form of kinetic energy, which results in a population of photoelectrons that have significantly larger energies than the thermal energy of the electron gas. These high energy photoelectrons then lose their energy by collisions with other atmospheric particles. For the thermal electrons, elastic collisions with non-thermal electrons is the main heat source and is the reason why the electron gas becomes hotter than the neutral and ion gases in the Earth’s upper thermosphere (Smithtro & Solomon 2008).

The two main assumptions in our model are that the photoelectrons lose their energy locally where they are created and that the non-thermal electron spectrum is in a steady state at each point. Given the latter assumption, the spectrum can be calculated simply by balancing sources and sinks of electrons at each electron energy. In future models, we will include also a more sophisticated electron transport model; this was shown to influence the heat deposition by Tian et al. (2008b). For the effects of collisions, the situation is complicated by the fact that a given atmospheric species can interact with a non-thermal electron in multiple ways. Each of these different interactions has a different energy-dependent cross-section and takes a different amount of energy from the impacting electron.

At a given electron energy, E_e, the two sources of electrons are photoionization reactions producing electrons with energy E_e, and the degradation of more energetic electrons through collisions with the ambient gas. The production spectrum for photoelectrons by photoionization reactions, P_e(E), is given by $$P_{\text{e}}(E_{\text{e}})={\displaystyle \sum _k\frac{I_{\text{xuv}}(E_{\text{e}}+\delta E_k)}{E_{\text{e}}+\delta E_k}}{\sigma }_k(E_{\text{e}}+\delta E_k)[R_k],$$(17)

where the sum is over all photoionization reactions considered, σ_k (E_e) is the energy dependent cross-section for the kth photoionization reaction, [R_k] is the number density of the reactant in the kth photoionization reaction, and δE_k is the energy required for the ionization to take place. We calculate the non-thermal electron spectrum using $${\varphi }_{\text{e}}(E_{\text{e}})=\frac{P_{\text{e}}(E_{\text{e}})+{\displaystyle \sum _i{\displaystyle \sum _j{n}_i}}{\sigma }_{ij}(E_{\text{e}}+\delta E_{ij}){\varphi }_{\text{e}}(E_{\text{e}}+\delta E_{ij})}{{\displaystyle \sum _i{\displaystyle \sum _j{n}_i}}{\sigma }_{ij}(E_{\text{e}})},$$(18)

where ϕ_e(E_e) is the electron flux at energy E_e and σ_ij is the cross-sections for electron impact interactions. This expression is described in more detail by Schunk & Nagy (1978). In both the numerator and the denominator, the first sum is over all species that the electrons interact with and the second sum is over all possible interactions with that species. The second term in the numerator is the source term from the degradation of higher energy electrons. The fact that photoelectrons can only lose energy, so that ϕ_e at a given energy depends on the higher energies values of ϕ_e only, makes solving Eq. (18) trivial. To do this, we break the spectrum down into 100 discreetenergy bins logarithmically spaced between 1 and 1000 eV. We first calculate ϕ_e at the bin with the highest energy assuming ϕ_e = 0 for higher energies, and then iterate downwards through the spectrum, calculating ϕ_e in each bin.

The neutral interacting species that we consider are N₂, O₂, O, CO₂, CO, and He. For N₂, we use the electron impact cross-sections given in Green & Barth (1965). For non-ionising O₂ transitions, we use the cross-sections from Watson et al. (1967), and for O₂ ionizations we use cross-sections from Jackman et al. (1977). For O, CO₂, and CO, we use cross-sections from Jackman et al. (1977). For He, we use cross-sections from Jusick et al. (1967).

In Fig. 4 we show our predicted non-thermal electron spectrum up to 70 eV for the current Earth’s atmosphere at an altitude of 200 km. This spectrum is an output of the Earth model presented in Sect. 3.1. For comparison, we also show the spectra calculated for several altitudes by Singhal & Haider (1984), which they compared to other models and observations in their Fig. 3. Our spectrum match theirs well for almost the entire energy range, indicating that our model calculates approximately realistic non-thermal electron spectra.

In our chemical network, we have included several ionization reactions due to impacts with non-thermal electrons. For each of these reactions, the total cross-sections for use in Eq. (21) are calculated by summing over the cross-sections for each corresponding ionization interaction. For many of the transitions that are not direct ionizations, an ionization can still take place by autoionization. We take these into account when calculating the total ionization cross-sections for O, CO₂, and CO by multiplying the cross-sections for these transitions by the autoionization factors given by Jackman et al. (1977). These reactions also remove energy from non-thermal photoelectrons, but the situation is complicated since the products of these reactions include two electrons. The energy of the original non-thermal electron that is not used to cause the reaction is distributed between the two electrons. For simplicity, we assume that one of the electrons gets this energy, and the other just becomes a normal thermal electron.

Fig. 4 Our prediction for the non-thermal electron spectrum at 200 km for the current Earth. For comparison, the dashed lines show the spectra at several altitudes calculated by Singhal & Haider (1984; taken from their Fig. 3 and multiplied by 4πto move the sr−1 from the units). The sudden increase in the flux at low electron energies is due to the thermal electron spectrum becoming dominant.

2.4 Chemical structure of the atmosphere

2.4.1 Chemistry

In this study, we attempt to construct a general chemical network that can be applied to a range of atmospheres with arbitrary compositions. This is difficult given the huge numbers of reactions and species that any network could consider and the uncertainties in the rate coefficients for the reactions, particularly at high temperatures. We do this by combining the networks of several previously published atmospheric models. The networks that we use are from Fox & Sung (2001), Verronen et al. (2002, 2005), Yelle (2004), García Muñoz (2007), Tian et al. (2008a), Richards & Voglozin (2011), Fox (2015), and Fox et al. (2015). We include almost all reactions from these papers, with some reactions being excluded if they introduced species that we consider unimportant. The rate coefficients are taken also from these studies in almost every case. Where multiple papers give different rate coefficients for the same reaction, we take the values almost arbitrarily, or find the coefficients on the KIDA database (Wakelam et al. 2012). In addition, we add a few reactions from KIDA that are not in any of these networks when necessary to stop reactions from creating species that are not destroyed. For the photoreactions, we take all of the relevant reactions from the PHIDRATES database (Huebner & Mukherjee 2015), which provides wavelength dependent cross-sections for the entire XUV spectrum. These cross-sections are all temperature independent, which is in many cases unrealistic and could lead to inaccuracies in our photochemistry (Venot et al. 2018). The few reactions involving non-thermal electrons produced in photoionization reactions are described in Sect. 2.3.3. The resulting network, which is given in Appendix H, contains 63 species, including 30 ion species, and 503 reactions, including 56 photoreactions and 7 photoelectron reactions.

The reaction rate of the kth chemical reaction, R_k, is related to the rate coefficient, k_k, by $$R_k=k_k{\displaystyle \prod _i{n}_i},$$(19)

where the RHS gives the product of the densities of all reactants. The rate coefficients for normal reactions are typically functions of temperature and many of the reactions have temperature limits, both of which are listed in Table H.1. A difficulty in our model is that we calculate separate neutral, ion, and electron temperatures, and in many cases it is unclear which of these temperatures to use to calculate the rate coefficients. For reactions that have only neutral reactants, we use the neutral temperature; for reactions that have a mixture of neutrals, ions, and electrons as reactants, we simply use the averages of the temperatures of the involved components (e.g. if a reaction has one neutral reactant and one ion reactant, we set T_gas = (T_n + T_i)∕2 in the equation for the rate coefficient). For photoreactions, the rate coefficients depend on the XUV spectrum and the wavelength dependent cross-sections by $$k_k={\displaystyle {\int }_{E_{\text{t}}}^{\infty }\frac{{\sigma }_k{I}_E}{E}}\text{d}E,$$(20)

where E is the photon energy, E_t is the threshold energy for the reaction, σ_k is the cross-section, and I_E is the irradience in units of energy flux per unitenergy (the quantity I_ν used elsewhere in this paper is the irradience in units of energy flux per unit frequency). Similarly, the equation for the rate coefficients of reactions involving inelastic collisions with non-thermal electrons is $$\begin{array}{lll}\hfill & \hfill & k_k={\displaystyle {\int }_{E_{\text{t}}}^{\infty }{\sigma }_k}{\varphi }_{\text{e}}\text{d}E.\hfill \end{array}$$(21)

The result is a set of ordinary differential equations, one for each species, describing the rates of change of the species densities. For the jth species, this can be written $$\frac{\text{d}{n}_j}{\text{d}t}={\displaystyle \sum _k}{R}_j-{\displaystyle \sum _i}{R}_i=S_j,$$(22)

where the first sum is over all reactions that create the jth species, and the second sum is over all reactions that destroy it. The S_j term is the total source term for the jth species in the RHS of Eq. (1). To evolve n_j using Eq. (22), we use an implicit Rosenbrock solver described in Appendix H.

In our model, we break the gas down into neutral, ion, and electron components. A difficulty in our model is that the chemical reactions cause the transfer of mass, momentum, and energy between the components simply due to the changes in the identities of atoms and molecules. For example, consider the reaction N⁺ + O₂ → NO⁺ + O; this reaction transfers an O atom from the neutral gas to the ion gas. The changes in the mass and momentum densities of the components are trivial to calculate, but the changes in the energy densities are not. We avoid this problem by assuming that the temperatures are unaffected when updating the species densities due to chemistry. The heating of the gas due to exothermic and endothermic chemical reactions, and the energy exchanges between the neutral, ion, and electron gases, are calculated separately, as described in Sect. 2.5.

2.4.2 Diffusion

Many of thespecies considered in the simulation are created and destroyed slowly by chemical and photochemical reactions. For these species, a very important transport mechanism is diffusion. Our model takes into account both molecular and eddy diffusion. Eddy diffusion evolves the density profiles so that they all follow the pressure scale height of the entire gas; molecular diffusion evolves the density profiles so that they all follow their own pressure scale heights. In the homosphere, eddy diffusion dominates and the mixing ratios of the long-lived species are independent of altitude. In the heterosphere, molecular diffusion dominates and the densities of heavy species decrease with increasing altitude faster than the densities of light species, meaning that light species become increasingly dominant at higher altitudes (this also happens due to the dissociation of heavy molecules). The equation that we use for the diffusive flux of the jth species, including both eddy and molecular diffusion, is $${\Phi }_{\text{d},j}=n_j{v}_{\text{d},j},$$(23)

where v_d,j is the diffusion speed, given by $$\begin{array}{l}\hfill v_{\text{d},j}=-D_j\left[\frac{1}{n_j}\frac{\text{d}{n}_j}{\text{d}r}-\frac{1}{N}\frac{\text{d}N}{\text{d}r}+\left(1-\frac{m_j}{\overline{m}}\right)\frac{1}{p}\frac{\text{d}p}{\text{d}r}+\frac{{\alpha }_{\text{T},j}}{T}\frac{\text{d}T}{\text{d}r}\right]\hfill \\ -K_{\text{E}}\left[\frac{1}{n_j}\frac{\text{d}{n}_j}{\text{d}r}-\frac{1}{N}\frac{\text{d}N}{\text{d}r}\right],\hfill \end{array}$$(24)

where D_j and K_E are the molecular and eddy diffusion coefficients, n_j and N are the particle number densities of the jth species and of the entire gas, m_j and $\overline{m}$ are the molecular masses of the jth species and of the entire gas, p and T are the thermal pressure and temperatureof the entire gas, and α_T,j is the thermal diffusion factor. To solve these equations, we use the implicit Crank–Nicolson method described in Appendix E.

For molecular diffusion, the diffusion coefficient for a given species, D_j, depends on both the species itself, the composition of the background gas, and the temperature. For all diffusion coefficients, we use the relation $$D_j=\frac{{\alpha }_j\times {10}^{17}{T}^{s_j}}{N}.$$(25)

For H, H₂, He, CH₄, CO, Ar, CO₂, and O, we use values for α_j and s_j given in Tables 15.1 and 15.2 of Banks & Kockarts (1973) assuming an N₂ background atmosphere, which is likely reasonable since the values are very similar for the other background atmospheres. For simplicity, we assume α_j = 1 and s_j = 0.75 for all other species. For H, H₂, and He, we assume α_T,j = −0.38 and for Ar, we assume α_T,j = 0.17 (Banks & Kockarts 1973); for all other species, we assume α_T,j = 0.

In models such as ours, the eddy diffusion coefficients as a function of altitude are free parameters; this is the only free parameter in our model. We assume it is given by K_E = AN^B, where K_E and N have the units cm² s⁻¹ and cm⁻³ respectively. This functional form is typically used for models of Venus and Mars (von Zahn et al. 1980; Fox & Sung 2001). For Venus, we use A = 2 × 10¹³ and B = −0.5 and impose a maximum value for K_E of 6 × 10⁸ cm² s⁻¹. These values were used by Fox (2015) for the upper atmosphere of Mars, and are very similar to values found for Venus (von Zahn et al. 1980). For the current Earth, we first fit A and B to the tabulated K_E values givenby Roble (1995), but we scale these values up by a factor of ten in order to fit the expected O₂ densities at high altitudes in our Earth model. This gives A = 10⁸ and B = −0.1. We note that even without scaling up the eddy diffusion coefficients from Roble (1995), we obtain good fits to the density profiles of all other species.

For atmospheres that are close to hydrostatic, it is important to specify a diffusion flux at the exobase. We assume an outward diffusion flux that corresponds to Jeans escape. This is only necessary for the lightest species (H and He), so for all species more massive than 4 m_p, we assume a zero flux. In simulations where the gas at the exobase is supersonic, and therefore is faster than the escape velocity, we simply assume a zero diffusion flux at the exobase.

2.5 Thermal structure of the atmosphere

2.5.1 Heating

In this model, the important energy sources are stellar XUV (10–4000 Å) and IR (1–20 μm) radiation. We also include a simple treatment of Joule heating. The total energy deposition rate from a radiation field travelling a distance dx through an absorbing gas is − dF∕dx, where F is the energy flux (=∫ I_νdν). In reality, the absorbed energy is not directly added to the thermal energy budget of the gas, and not all of the energy deposited is eventually converted to heat. A common way to calculate the heating is to multiply the total energy deposition rate by a heating efficiency factor (e.g. Erkaev et al. 2013; Johnstone et al. 2015b); this assumption is generally undesirable since it adds an unconstrained free parameter into the model. We instead use a more complete heating model where the energy release from different processes are calculated individually. The heating processes considered are direct heating by the stellar XUV field, electron heating by elastic collisions with non-thermal electrons, heating from exothermic chemical reactions, direct heating by the stellar IR field, and Joule heating. A diagram demonstrated how stellar XUV energy in our model is transferred to heat after been absorbed is shown in Fig. 5.

When a XUV photon is absorbed, a large part of its energy is used to cause either a dissociation, an ionization, or both. For photodissociation reactions, the remaining photon energy is given to the products as kinetic energy, and ultimately dissipated as heat in the gas. It is this heating that we consider the direct heating by the stellar XUV field. To calculate this direct heating, we need to consider each photodissociation reaction and each energy bin in the XUV field separately. The full equation for the heating rateis $$Q_{\text{xuv}}={\displaystyle \sum _k{Q}_{\text{xuv},\text{k}}},$$(26)

where Q_xuv,k is the heating rate by the kth reaction and is given by $$Q_{\text{xuv},\text{k}}={\displaystyle {\int }_{E_{\text{T},k}}^{\infty }\left(E-E_{\text{T},k}\right)}\frac{I_E}{E}{\sigma }_k(E)[R_k]\text{d}E,$$(27)

where the sum is over all photodissociation reactions, E_T,k is the energy required for the reaction to take place, σ_k(E) is the reaction cross-section at energy E, [R_k] is the number density of the reactant, and I_E is the irradience in units of energy flux per unit energy. The term I_E ∕E is the photon flux per unit photon energy at energy E and the term (I_E∕E)σ_k(E)[R_k] is the rate at which the kth reaction takes place per unit volume per unit photon energy. The energy released per reaction is E− E_T,k, which when multipled by the reaction rate and integrated over all photon energies gives the heating rate for the kth reaction.

For photoionization reactions, we do not give the remaining photon energy to the gas as heat directly, but instead assume that all this energy is given to the resulting free electron and we calculate the non-thermal electron spectrum, as described in Sect. 2.3.3. This energy is either given to the neutral gas by inelastic collisions, typically exciting atoms/molecules or causing secondary ionizations, or it is given to the thermal electrons by elastic collisions. We assume that the energy given to the neutral gas is all lost by radiative relaxation and consider therefore only the heating of the electron gas. Using the expression given by Schunk & Nagy (1978), we calculate the electron heating rate as $$\begin{array}{lll}{Q}_{\text{e}}\hfill & =\hfill & {\displaystyle {\int }_0^{E_{\text{t}}}\left(E_{\text{e}}-\frac{3}{2}{k}_{\text{B}}{T}_{\text{e}}\right)}{P}_{\text{e}}(E_{\text{e}})\text{d}{E}_{\text{e}}+{\displaystyle {\int }_{E_{\text{t}}}^{\infty }{n}_{\text{e}}}{L}_{\text{e}}(E_{\text{e}}){\varphi }_{\text{e}}(E_{\text{e}})\text{d}{E}_{\text{e}}\hfill \\ \hfill & \hfill & +\left(E_{\text{t}}-\frac{3}{2}{k}_{\text{B}}{T}_{\text{e}}\right)n_{\text{e}}{L}_{\text{e}}(E_{\text{t}}){\varphi }_{\text{e}}\left(E_{\text{t}}\right),\hfill \end{array}$$(28)

where E_e is the electron kinetic energy, P_e(E_e) is the production spectrum of electrons (see Eq. (17)), E_t is the energy above which the non-thermal flux is larger than the thermal flux, and L_e (E_e) is the loss function given by $$L_{\text{e}}(E_{\text{e}})=\frac{3.37\times {10}^{-12}}{E_{\text{e}}^{0.94}{n}_{\text{e}}^{0.03}}{\left(\frac{E_{\text{e}}-E_{\text{th}}}{E_{\text{e}}-0.53E_{\text{th}}}\right)}^{2.36},$$(29)

where E_th = 8.618 × 10⁻⁵ T_e (Swartz et al. 1971). The three terms on the RHS of Eq. (28) give respectively the heating/cooling by the direct production of thermal electrons by photoionization reactions, heating of thermal electrons by elastic collisions with non-thermal electrons, and a surface term related to the crossover between the thermal and non-thermal electron spectra. Although a more accurate version of the surface term was derived by Hoegy (1984), we use the version given above because it is simpler and is sufficiently accurate for our purposes.

Much of the photon energy used to cause a photoreaction is not lost, but is instead converted into chemical potential energy that can then be released as heat in exothermic chemical reactions. The heating rate at a given point by chemical reactions is given by $$Q_{\text{chem}}={\displaystyle \sum _k{R}_k}{Q}_{\text{chem},k},$$(30)

where the sum is over all chemical reactions, and R_k and Q_chem,k are the reaction rate and energy released per reaction for the kth reaction. The values of Q_chem,k are given for each reaction in Table H.1. These energies are mostly taken from Tian et al. (2008a) or from the KIDA database, and when the energy for a given reaction is not available in either of these sources, we simply assume it does not contribute to the heating.

We consider also heating of the atmosphere by the abosrption of IR radiation. We assume that all of the energy removed from the IR spectrum is input into the neutral gas as thermal energy, giving a heating rate of $$Q_{\text{IR}}={\displaystyle \sum _j{\displaystyle {\int }_{\nu }{\sigma }_{\nu ,j}}}[R_j]I_{\nu }\text{d}\nu ,$$(31)

where the sum is over all absorbing species, the integral is over the entire IR spectrum that we consider, and σ_ν,j and [R_j] are the absorption cross-section and number density of the jth absorbing species. In reality, this energy is first used to excite CO₂ and is then released as heat through collisional deexcitation, which we take into account with an additional excitation term in the equations for CO₂ cooling.

Additionally, the upper atmospheres of magnetised planets are heated by two magnetospheric processes: these are energetic particle precipitation and Joule heating. For the Earth, during quiet geomagnetic conditions these two processes are likely similar in magnitude, and Joule heating tends to dominate during geomagnetic storms (e.g. Chappell 2016). This process could become important for planets that are exposed to extreme space weather (Cohen et al. 2014). Both processes are most significant at high latitudes, but tend not to influence the global heat budget significantly during quiet conditions. In this paper, we model Joule heating using the simplified model described in Roble et al. (1987) and Smithtro & Sojka (2005). The two input parameters are the ambient magnetic field strength, which we assume is 0.5 G everywhere, and the total global Joule heating rate, which we assume is 1.4 × 10¹⁸ erg s⁻¹. This is double the value used in Roble et al. (1987), which is typical for quiet levels of geomagnetic activity (Foster et al. 1983); we double the value to take into account also the energy input expected from particle precipitation. We calculate the heating rate at each altitude using $$Q_{\text{J}}={\sigma }_{\text{P}}{E}^2,$$(32)

where E is the electric field strength and σ_P is the Pedersen conductivity (Foster et al. 1983). We do not calculate the electric field, but instead assume that it is a constant and use it as a free parameter that can be scaled in order to give us the desired total global Joule heating rate. The Pedersen conductivity varies with altitude, and at a given point depends on the densities of individual ion and neutral species, the gas temperature, and the ambient magnetic field strength. The σ_P profiles are calculated self-consistently within the model using the equations described in Sect. 5.11 of Schunk & Nagy (2000). The equation for σ_P is $${\sigma }_{\text{P}}={\displaystyle \sum _i{\sigma }_i}\frac{{\nu }_i^2}{{\nu }_i^2+{\omega }_i^2},$$(33)

where the sum is over all considered ion species, and σ_i, ν_i, and ω_i are the ion conductivity, ion-neutral collision frequency, and angular gyrofrequency of the ith ion species. The angular gyrofrequency is given by ω_i = q_iB∕m_i. The ion conductivity is given by ${\sigma }_i=(n_i{q}_i^2)/(m_i{\nu }_i)$, where n_i, q_i, and m_i are the ion number density, charge, and mass respectively. The ion-neutral collision frequency, ν_i, is calculated as the sum over the collision frequencies with individual neutral species, such that ${\nu }_i={\displaystyle \sum _n{\nu }_{in}}$, where ν_in is the frequency of collisions between the ith ion species and the nth neutral species. For this, we use the same collisions and ν_in values described in Sect. 2.5.4 for ion-neutral heat exchange.

Fig. 5 Simplified cartoon illustrating the main pathways taken in our model by the energy that is removed from the XUV field by absorption due to photodissociation and photoionization reactions. In both cases, much of the energy used to cause the photoreaction is released due to exothermic chemical reactions. For photodissociation, the remaining energy from the absorbed photon is given to the thermal energy budget of the neutral gas directly. For photoionization, the remaining energy is released as kinetic energy of the produced electron, and is then lost as the electron collides with the ambient gas; collisions with neutral species are inelastic and lead to excitation, dissociation, and ionization, whereas collisions with ambient thermal electrons are elastic and lead to heating of the electron gas.

2.5.2 Cooling

We consider the effects of IR cooling by CO₂, H₂ O, NO, and O. In all cases, cooling happens when atoms/molecules are excited by collisions with other particles and then radiate the energy awaybefore they are deexcited by further collisions. Collisions cause there to be a continuous transfer of energy from the atmosphere’s thermal energy reservoir to the various forms of energy within the individual atoms and molecules, and a corresponding transfer of energy back to the thermal energy reservoir. However, due to radiative relaxation (i.e. spontaneous/stimulated emission) and the loss of many of the emitted photons to space, the rate at which energy is transferred back to the thermal reservoir is reduced, and the resulting inbalance is the cooling that we are interested in. The calculation of the cooling rates is seldom trivial and ideally would involve calculating the full transport ofthe emitted IR spectrum through the atmosphere and tracking the populations of each of the various excited states in the relevant species (e.g. see Wintersteiner et al. 1992). In this paper, we take into account all of these processes in a simpler way and aim to implement more sophisticated treatments of cooling in future studies.

Cooling by CO₂ is dominated by emission at 15 μm. We use the cool-to-space approximation (e.g. Dickinson 1972); the fundamental assumption is that the cooling at each point is caused entirely by emitted photons that escape directly to space. Ignoring stimulated emission, this means $$Q_{{\text{CO}}_2}={(h\nu )}_{15\mu \text{m}}{A}_{10}[{\text{CO}}_2^{*}]ϵ,$$(34)

where ${(h\nu )}_{15\mu \text{m}}$ is the energy of a single 15 μm photon (=1.325 × 10⁻¹³ erg), A₁₀ is the Einstein coefficient for spontaneous emission, $[{\text{CO}}_2^{*}]$ is the density of excited CO₂ molecules, and ϵ is the probability that a 15 μm photon emitted from a given point escapes to space. The definition of ϵ is such that it takes into account the fact that photons emitted downwards are not lost, and therefore approaches a maximum of 0.5 at high altitudes.

In order to calculate $[{\text{CO}}_2^{*}]$, we consider three excitation and two deexcitation mechanisms. The excitation mechanisms are collisional excitation, the absorption of 15 μm photons previously emitted by excited CO₂ molecules, and the absorption of photons from the host star’s IR spectrum. For the second process, the fundamental simplifying assumption is that all photons that are emitted and are not lost to space are reabsorbed locally where the emission took place. The excitation rate is given by $$\begin{array}{lll}\frac{\text{d}[{\text{CO}}_2^{*}]}{\text{d}t}\hfill & =\hfill & \left({\displaystyle \sum _M{k}_{e,M}}[M]\right)\left([{\text{CO}}_2]-[{\text{CO}}_2^{*}]\right)\hfill \\ \hfill & \hfill & +A_{10}[{\text{CO}}_2^{*}]\left(1-ϵ\right)+S_{\text{IR}},\hfill \end{array}$$(35)

where the sum is over all species that collisionally excite CO₂, k_e,M is the rate coefficient for collisional excitation, and S_IR is the additional excitation term due to the absorption of stellar IR radiation given by Eq. (16). The term $[{\text{CO}}_2]-[{\text{CO}}_2^{*}]$ is the density of non-excited CO₂ molecules. The deexcitation mechanisms are collisional deexcitation and radiative relaxation. The deexcitation rate is given by $$-\frac{\text{d}[{\text{CO}}_2^{*}]}{\text{d}t}=\left({\displaystyle \sum _M{k}_{d,M}}[M]\right)[{\text{CO}}_2^{*}]+A_{10}[{\text{CO}}_2^{*}].$$(36)

Where the first term on the RHS dominates, the atmosphere is in the local thermodynamic equilibrium (LTE) regime, and where the two terms are similar, or the second term dominates, the atmosphere is in the non-local thermodynamic equilibrium (non-LTE) regime. Assuming a steady state, Eqs. (35) and (36) add up to zero, giving $$\left[{\text{CO}}_2^{*}\right]=\frac{{\displaystyle \sum _M{k}_{e,M}}[M][{\text{CO}}_2]+S_{\text{IR}}}{{\displaystyle \sum _M\left(k_{e,M}+k_{d,M}\right)}[M]+A_{10}ϵ}.$$(37)

The Einstein coefficient, A₁₀, is 0.46 s⁻¹ (Curtis & Goody 1956), and the rate coefficients are related by $k_{e,M}=2k_{d,M}\exp \left(-667/T_{\text{n}}\right)$ (Castle et al. 2006). For the escape probabilities, we use the tabulated values given by Kumer & James (1974) which depend entirely on the amount of CO₂ above the considered altitude, z, given by $N_{\text{c}{\text{o}}_{\text{2}}}={\int }_z^{\infty }\left[\text{C}{0}_2\right]\text{d}z$. We fit their tabulated values with $$ϵ=\{\begin{array}{ll}0.7202{\left(\sigma N_{{\text{CO}}_{\text{2}}}\right)}^{-0.613}\hfill & \text{if}\sigma N_{{\text{CO}}_{\text{2}}}>2,\hfill \\ 0.4732{\left(\sigma N_{{\text{CO}}_{\text{2}}}\right)}^{-0.0069}\hfill & \text{if}\sigma N_{{\text{CO}}_{\text{2}}}<2.\hfill \end{array}$$(38)

where σ =6.43 × 10⁻¹⁵ cm². The dependence of ϵ on $\sigma N_{{\text{CO}}_{\text{2}}}$ is shown in Fig. 6. For the collisional excitation/deexcitation rate, we consider the influences of O, O₂, N₂, CO₂, He, and Ar. For O, we use the experimentally measured values of k_d,M given by Castle et al. (2012), and for the other species, we use the measured values given by Siddles et al. (1994). In all cases, k_d,M have temperature dependences that we fit using power-laws of the form $k_{d,M}(T_{\text{n}})=A{T}_{\text{n}}^B$, where the values of A and B are given in Table 1. In Fig. 6, we show the measured deexcitation rates and our analytic fit formulae for each species. A signficant worry with our fit formulae for k_d,M is that all of the measurements that we use are for low gas temperatures, and therefore our fit formulae might be inaccurate at high temperatures. This problem is not likely to influence our results in this paper since CO₂ cooling is only significant in regions of the atmospheres that are within the experimental temperature ranges.

For NO cooling, we use the model given by Oberheide et al. (2013). We consider emission in the vibrational band at 5.3 μm assuming two excitation mechanisms: these are collisional excitation by O atoms and radiative pumping by earthshine. We assume that all photons emitted by NO molecules escape to space (i.e. ϵ = 1), which is realistic for the Earth since NO cooling is only significant in the thermosphere (Kockarts 1980). The cooling rate is given by $$Q_{\text{NO}}={(h\nu )}_{5.3\mu \text{m}}{A}_{10}[{\text{NO}}^{*}],$$(39)

where ${(h\nu )}_{5.3\mu \text{m}}=3.75\times {10}^{-13}$ erg and [NO^*] is the density of excited NO molecules, given by $$[{\text{NO}}^{*}]=\frac{k_{e,\text{O}}[\text{O}]+S_{\text{E}}}{\left(k_{e,\text{O}}+k_{d,\text{O}}\right)[\text{O}]+S_{\text{E}}+A}[\text{NO}],$$(40)

where S_E is the excitation rate due to earthshine, k_e,O and k_d,O are the collisional excitation and deexcitation rate coefficients, and A is the Einstein coefficient for spontaneous emission. As in Oberheide et al. (2013), we use S_E = 1.06 × 10⁻⁴ s⁻¹, k_d,O = 2.8 × 10⁻¹¹ cm³ s⁻¹, A =12.54 s⁻¹, and $k_{e,\text{O}}=k_{d,\text{O}}\exp \left(-2700/T_{\text{n}}\right)$.

For O cooling, we consider emission at 63 μm and 147 μm using the parameterization derived by Bates (1951; see Eqs. (14.57) and (14.58) of Banks & Kockarts 1973) given by $$Q_{\text{O}}=Q_{\text{O},63\mu \text{m}}+Q_{\text{O},147\mu \text{m}},$$(41)

where $$\begin{array}{lll}\hfill & \hfill & Q_{\text{O},63\mu \text{m}}=\frac{1.67\times {10}^{-18}\exp \left(-228/T_{\text{n}}\right)[\text{O}]}{1+0.6\exp \left(-228/T_{\text{n}}\right)+0.2\exp \left(-326/T_{\text{n}}\right)},\hfill \\ \hfill & \hfill & Q_{\text{O},147\mu \text{m}}=\frac{4.59\times {10}^{-20}\exp \left(-326/T_{\text{n}}\right)[\text{O}]}{1+0.6\exp \left(-228/T_{\text{n}}\right)+0.2\exp \left(-326/T_{\text{n}}\right)},\hfill \end{array}$$

where [O] is in cm⁻³, T_n is in K, and the cooling rates are in erg s⁻¹ cm⁻³. More sophisticated modelling of O cooling will be used in future models.

For cooling by H₂O, we use the parametitzation for emission in rotational bands by Hollenbach & McKee (1979) and summarised in Kasting & Pollack (1983, see their Eqs. (32)–(38)). In this model, H₂ O is excited by collisions with H atoms only. Given the length of the set of equations involved, we do not write them here.

Fig. 6 Escape probability of a 15 μm photon as a function of CO2 optical depth (upper panel) and the CO2 deexcitation rate coefficients as a function of temperature (lower panel). In the upper panel, the numbers show the approximate altitudes in our current Earth model where these points on the line are reached. In the lower panel, the data points are measurements from Siddles et al. (1994) and Castle et al. (2012), and the solid lines are our power-law fits, given in Table 1.

2.5.3 Conduction

In the Earth’s upper thermosphere, cooling of the neutral gas is not strong enough to balance heating, and a steady state is only reached because conduction downwards into the cooler lower thermosphere removes this excess energy. Since the temperatures of the neutrals, ions, and electrons are evolved separately, separate conductivities must be used for each of these components. For ion and electronconductivities, we ignore the effects of the magnetic field, which reduces the conduction in directions perpendicular to the magnetic field. The conduction equations are solved using the implicit Crank–Nicolson method, as described in Appendix F.

For the neutral gas, we consider eddy conduction, which is dominant in the lower atmosphere, and molecular conduction, which is dominant in the upper atmosphere. The neutral conduction equation is $$\frac{\partial e_{\text{n}}}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2{\kappa }_{\text{mol}}\frac{\partial T_{\text{n}}}{\partial r}+r^2{\kappa }_{\text{eddy}}\left(\frac{\partial T_{\text{n}}}{\partial r}+\frac{g}{c_{\text{P}}}\right)\right],$$(44)

where κ_mol is the molecular conductivity, κ_eddy is the eddy conductivity, g is the gravitational acceleration, and c_P is the specific heat at constant pressure. The term g∕c_P for the eddyconduction is the adiabatic lapse rate. The eddy conductivity is related to the eddy diffusion coefficient by κ_eddy = ρc_PK_E (Hunten 1974). The molecular conductivity is dependent on the temperature and composition of the gas. We estimate κ_mol using the equations given in Sect. 14.3 of Banks & Kockarts (1973) with some minor simplifications. The molecular conductivity of the kth species is given by $${\kappa }_k=A_k{T}_{\text{n}}^{s_k},$$(45)

where A_k and s_k are coefficients that depend on the species. We assume the total conductivity of the gas is given by $${\kappa }_{\text{mol}}={\displaystyle \sum _k\frac{n_k{\kappa }_k}{{\displaystyle \sum _j{n}_j}{\varphi }_{kj}}},$$(46)

where $${\varphi }_{kj}=\frac{{\left[1+{\left({\kappa }_k/{\kappa }_j\right)}^{\frac{1}{2}}{\left(m_j/m_k\right)}^{\frac{1}{4}}\right]}^2}{2\sqrt{2}{\left[1+\left(m_j/m_k\right)\right]}^{\frac{1}{2}}},$$(47)

where m_k is the molecular mass of the kth species. Thesums in Eq. (46) should be over all neutral species, but in reality we only consider species for which we have A_k and s_k values. The species we consider are N₂, O₂, CO₂, CO, O, He, H, and Ar, with values for A_k and s_k taken from Table 13 of Bauer & Lammer (2004) for Ar, and Table 10.1 of Schunk & Nagy (2000) for the others.

For the ions, the conduction equation is $$\frac{\partial e_{\text{i}}}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2{\kappa }_{\text{i}}\frac{\partial T_{\text{i}}}{\partial r}\right],$$(48)

where κ_i is the ion conductivity. For κ_i, we use Eq. (22.122) from Banks & Kockarts (1973) which expresses the conductivity of ion gases made of a single ion species as $4.6\times {10}^4{A}^{-1/2}{T}_{\text{i}}^{5/2}$ eV cm⁻¹ s⁻¹ K⁻¹, where A is the atomic mass of the species and T_i is the ion temperature. For a gas mixture, they recommend using a density weighted average thermal conductivity, so we adopt the form $${\kappa }_{\text{i}}=4.6\times {10}^4\frac{{\sum }_k{n}_k{A}_k^{-\frac{1}{2}}}{{\sum }_k{n}_k}{T}_{\text{i}}^{\frac{5}{2}},$$(49)

where the sums are over all ion species and n_k is in cm⁻³, T_i is the ion temperature in K, and κ_i is in eV cm⁻¹ s⁻¹ K⁻¹.

For the thermal electron gas, the conduction equation is the same as Eq. (48) with the subscript i replaced by e. We calculate the electron conductivity using Eq. (22.116) of Banks & Kockarts (1973): $${\kappa }_{\text{e}}=\frac{7.7\times {10}^5{T}_{\text{e}}^{5/2}}{1+3.22\times {10}^4\left(T_{\text{e}}^2/n_{\text{e}}\right){\displaystyle \sum _k{n}_k}{\overline{Q}}_{\text{D},k}},$$(50)

where ${\overline{Q}}_{\text{D},k}$ is the average momentum transfer cross-section of the kth species. The sum in the denominator should technically be over all neutral species, but in reality only the main species contribute significantly. In this sum, we take into account the effects of N₂, O₂, O, H, and He using the temperaturedependent equations for ${\overline{Q}}_{\text{D},k}$ given in Table 9.2 of Banks & Kockarts (1973). The numerator in Eq. (50) is the electron conductivity of a fully ionised gas, and the denominator corrects for the reduction in conductivity caused by collisions with neutrals reducing the mean free paths of thermal electrons.

2.5.4 Energy exchange

The neutral, ion, and electron gases exchange energy by collisions. In our model, the electrons lose energy only by collisions with neutrals and ions. The energy gained by the ions is then given to the neutrals by further collisions, which is the most important neutral heating mechanism in the upper thermosphere. The energy exchange equations are solved using the implicit Crank–Nicolson method, as described in Appendix G.

For the electron-ion energy exchange, we take into account elastic Coulomb collisions only and the total energy exchange rate is calculated by summing over the rates for individual ion species. The basic equation is $$Q_{\text{ei}}=-3k_{\text{B}}\left(T_{\text{e}}-T_{\text{i}}\right){\sum }_k\frac{n_k{m}_k{\nu }_{\text{e}k}}{m_{\text{e}}+m_k},$$(51)

where the sum is over all ion species and ν_ek is the momentum transfer collision frequency between electrons and the kth ion species, This equation, derived by Schunk (1975), requires several assumptions, including that the temperature difference between the electrons and ions are small. The definition of Q_ei is such that a positive value means energy is taken from the ions and given to the electrons. To calculate ν_{e k} for a given ion, we use ν_{e k} = 54.5 $n_k{Z}_k^2/T_{\text{e}}^{3/2}$, where Z_k is the charge of the ion (see Sect. 4.8 of Schunk & Nagy 2000).

The total ion-neutral energy exchange rate is the sum of the exchange rates of individual species pairs. The equation for the energy exchange rate is $$Q_{\text{in}}=-3k_{\text{B}}\left(T_{\text{i}}-T_{\text{n}}\right){\sum }_n{\sum }_k\frac{n_i{m}_i{\nu }_{in}}{m_i+m_n},$$(52)

where the subscripts n and i are for the nth neutral and ith ion species. The definition of Q_in is such that a positive value means energy is taken from the neutrals and given to the ions. The neutrals that we consider in the sums are H, He, N, O, CO, N₂, O₂, and CO₂; the ions that we consider are H⁺, He⁺, C⁺, N⁺, O⁺, CO⁺, N${}_2^{+}$, NO⁺, O${}_2^{+}$, and CO${}_2^{+}$. The ion-neutral heat exchange, described in detail in Schunk & Nagy (2000), is dominated by two types of interactions: resonant and non-resonant interactions which dominate at high (>300 K) and low temperatures respectively. The resonant interactions happen when a neutral approaches its ion equivalent (e.g. O and O⁺) and charge exchanges with it, with the changing identities of the particles representing a net energy exchange between the ion and neutral gases. We use the temperature-dependent equations for the momentum transfer collision frequencies, ν_in, for individual ion-neutral pairs given in Table 4.5 of Schunk & Nagy (2000). Non-resonant interactions involve neutral species and dissimilar ions. For these interactions, ν_in can be described simply by ν_in = C_inn_n, where we use the coefficients C_in for individual ion-neutral pairs given in Table 4.4 of Schunk & Nagy (2000).

Several important mechanisms exist that cause thermal electrons to lose energy to the neutral gas. Electron–neutral heat exchange is very important for the electron temperature structure in the low thermosphere of the Earth, and normally proceeds through inelastic collisions that excite neutral atoms or molecules. The most important of these processes, at least for the current Earth, is inelastic collisions that cause fine structure transitions in ground state atomic oxygen; for this process, we use the scaling laws derived by Hoegy (1976). We also consider the excitation of ground state oxygen to the O(¹ D) excited state. In addition, we take into account energy exchange from electron collisions with neutral molecules that cause the exitation of rotational or vibrational modes in the molecule; for these, we consider collisions with N₂, O₂, H₂, CO₂, CO, and H₂O. We use the scaling laws for these processes that are convieniently listed in Sect. 9.7 of Schunk & Nagy (2000).

3 Model validation

To validate our model, we calculate the upper atmospheres of modern Earth and Venus in this section. In both simulations, we use the modern solar XUV spectrum given by Claire et al. (2012), which represents the Sun approximately at the maximum of its activity cycle.

3.1 Earth

To validate our model for the Earth, we compare our results to those of two empirical models. For the neutral gas, we use the atmospheric density and temperature profiles of the empirical NRLMSISE-00 model (Picone et al. 2002). This model produces vertical profiles for the Earth’s atmosphere at arbitrary longitudes and latitudes and at arbitrary dates; the output profiles that we use are for temperature, and the densities of N₂, O₂, N, O, H, Ar, and He. For the ion densities, we use International Reference Ionosphere 2007 (IRI-2007; Bilitza & Reinisch 2008). This model produces vertical profiles for O${}_2^{+}$, NO⁺, O⁺, N⁺, H⁺, and electrons. We use these two standard models to obtain vertical atmospheric profiles for all longitudes and latitudes on the 1st January 1990, when the Sun was at approximately peak activity. We then calculate globally averaged profiles for this date. For our own Earth simulation, we assume a zenith angle of 66°, which we show below provides a good approximation for the globally averaged profiles.

We model the Earth’s upper atmosphere between an altitude of 65 km and the exobase. At the lower boundary, we use the values for temperature and density from this altitude in the NRLMSISE-00 model for comparison purposes, and additionally assume CO₂ and H₂ O mixing ratios of 4 × 10⁻⁴ and 6 × 10⁻⁶ respectively, which are reasonable values for the Earth’s middle atmosphere (Körner & Sonnemann 2001). In Fig. 7, we show the thermal structure of our current Earth model. The dashed line shows the standard atmosphere model that we use for comparison for the neutral gas temperature. The comparison between our results and the standard atmosphere model is very good, though an exact match between the models should not be expected especially since our input solar XUV spectrum will not match exactly the one used to produce the standard model. In Fig. 8, we show the strengths of the various heating and cooling mechanisms for this model. The results resemble very closely those of other global upper atmosphere models (e.g. Roble 1995; Tian et al. 2008a).

In Fig. 9, we show the densities as a function of altitude of several important species in our simulation. The species are chosen to be those output by the NRLMSISE-00 and IRI-2007 models that we use for comparison. Our predicted density structures are very similar to those of the comparison models, with the only exception being He, which we predict to be less abundant at high altitudes than expected. Clearly our model is able to realistically predict the structure of the Earth’s atmosphere.

Fig. 7 Neutral, ion, and electron temperature structures of our current Earth model. The dashed black line is for the empirical NRLMSISE-00 model.

Fig. 8 Thermal processes in our Earth model as functions of altitude. In the upper panel, the red line shows the total XUV, chemical, and Joule heating, the blue line shows the total IR cooling, the cyan and magenta lines show the neutral heating by collisions with electron and ions, the black line shows the sum of all of these, and the green line shows the effects conduction. In the middle and lower panels, we show the contributions of individual mechanisms to the total heating and cooling.

Fig. 9 Density structures of several important neutral (upper panel) and ion (lower panel) species. The solid lines show the predictions of our model and the dotted lines are from the empirical NRLMSISE-00 model for the neutrals and the IRI-2007 model for the ions.

3.2 Venus

In this section, we further validate our model by calculating the structure of the upper atmosphere of Venus. This is especially useful since Venus’ atmosphere is made up mostly of CO₂ and has therefore much stronger atmospheric cooling, allowing us to test that our model realistically responds to large changes in the CO₂ content. For a reference atmosphere, we use the empirical thermosphere model given in Table 3b of Hedin et al. (1983), which is for Venus’ atmosphere at noon during approximately solar maximum conditions. We therefore assume a zenith angle of 0° in our model and use the values given by Hedin et al. (1983) for the temperature and species densities at the lower boundary. We do not include Joule heating in our Venus model.

In Fig. 10, we show a summary of the results of our Venus model. We compare our calculated temperature structures to the standard model given by Hedin et al. (1983) and to the recent 3D global models for the full atmosphere of Venus by Gilli et al. (2017). The neutral temperature profile resembles that of Hedin et al. (1983), with very similar exobase temperatures, though we find that our model is colder in the lower thermosphere, with the largest difference being around 40 K around 130 km. These differences in the temperature profiles do not suggest that there is a problem with our simulations; the Gilli et al. (2017) model is also colder than the Hedin et al. (1983) profiles at 140 km, though it is much warmer than both our model and the Hedin et al. (1983) lower in the thermosphere. In Fig. 10, the density profiles for several species from our simulations can be compared to those from Hedin et al. (1983). As with the case of the Earth’s thermosphere, we underestimate the He abundances at high altitudes. Our other density profiles are similar to those of Hedin et al. (1983), with the differences being mostly due to the different temperatures.

Fig. 10 Results of our modern Venus model. In the upper left panel, showing the temperature profiles, the dashed black line shows the empirical model by Hedin et al. (1983) and the dotted black line shows the daytime profiles shown in Fig. 9 of Gilli et al. (2017). In the upper right panel, showing the density profiles of several important species, the solid lines show the results of our model and the dotted lines show the profiles given by Hedin et al. (1983). The lower left panel shows the heating and cooling mechanisms for the neutral gas, and the lower right panel shows all of the heating mechanisms.

4 Results

We present two applications of our code. In Sect. 4.1, we explore the effects of enhancing the CO₂ abundance in the current Earth’s atmosphere on the upper atmospheric structure. In Sect. 4.2, we explore the response of Earth’s upper atmosphere to the evolving XUV spectrum of the Sun.

4.1 The influence of enhanced CO₂ abundance on the atmospheric structure

To study the response of the Earth’s upper atmosphere to large changes in the CO₂ abundance, we calculate several models for the current Earth, varying only the lower boundary density of CO₂. We calculate models where the CO₂ base density is 0.01, 0.1, 1, 10, 100, and 1000 times the value for the current Earth. We refer to these models by this multiplicative factor, such that our model with 1000 times more CO₂ is the 1000× model. In our 1000× model, CO₂ has a mixing ratio of approximately 0.3, which is a factor of a few lower than that of Venus. As in previous sections, our input solar spectrum is for the current Sun at activity maximum.

It is important to note that the fact that we do not include the effects of changing the CO₂ abundances on the lower atmosphere. In reality, large changes in the CO₂ abundances will cause also changes in these lower boundary properties. This has the effect that the total heating and cooling rates approximately balance only in our 1× model. In models with lower CO₂ abundances, energy is removed from the computational domain at the lower boundary by downward conduction, causing the total cooling to be less than the total heating. In models with higher CO₂ abundances, the opposite effect takes place.

In Fig. 11, we show the neutral temperature structures for each model. As expected, increasing the CO₂ abundance leads to a significant decrease in the thermospheric temperature due to enhanced CO₂ cooling. Similarly, decreasing the amount of CO₂ leads to the thermosphere becoming hotter. In the 1000 × model, the maximum temperature reached in the thermosphere is only ~300 K, similar to that of Venus, and the mesosphere is cooled to a minimum temperature of ~130 K, which is much cooler than the minimum temperature of ~200 K that we find in our current Earth model. Fig. 11 also shows the total atmospheric heating rate⁴ as a function of CO₂ mixing ratio for different heating processes. We note however that these quantities are relatively crude measures of how important the various heating mechanisms are; for example, the total photoelectron heating for the current Earth is very small, but the effect is relatively large because it takes place high in the thermosphere where the gas densities are low and therefore less energy is required for the heating to be significant.

In the current Earth case, the two dominant heating mechanisms are chemical and XUV heating. These two mechanisms depend primarily on the input XUV energy flux, and therefore remain approximately constant. Due to the decrease in the ionization fraction of the gas, the heating of the electrons by collisions with non-thermal photoelectrons decreases significantly as the CO₂ mixing ratio is increased (see the linear dependence between electron heating and density in Eq. (28)). The mechanism that changes the most is the heating due to the absorption of stellar IR photons. Unlike the stellar XUV photons, which are all absorbed in the upper atmosphere in all cases, most of the IR photons pass through the upper atmosphere unhindered; therefore, adding efficient IR absorbing gases influences significantly how much IR energy is absorbed. In the 1000× model, we find that most of the heating is from IR absorption, especially at low altitudes. In the 0.01× model, the IR heating is dominated by H₂O absorption.

It is interesting to compare our simulations to those of Kulikov et al. (2007) who also modelled the effects on enhanced CO₂ abundances on the Earth’s thermosphere (see their Fig. 2). Our results are broadly similar which provides additional validation of our model. We find in general cooler upper thermospheric temperatures than they do for the enhanced CO₂ models. Our models differ in a few important ways: for example, they did not calculate separate neutral, ion, and electron temperatures and they did not consider the effects of chemical reactions on the density structures of individual species. The differences between our results likely have two sources. Firstly, the XUV heating in their model is based on an assumed heating efficiency parameter that in reality might vary with CO₂ abundance. Secondly, they put the lower boundaries of their models at the mesopause at an altitude of 100 km with an approximately fixed temperature, whereas we calculate the upper mesosphere from 65 km. In the enhanced CO₂ models, we get significant additional mesospheric cooling, and therefore much lower mesopause temperatures, which leads to cooler temperatures also at higher altitudes.

The changes in the exobase are demonstrated in Fig. 12. In the models with larger CO₂ abundances, the exobase altitudes and gas temperatures are much lower. Interestingly, the electron temperature does not decrease as much as the neutral and ion temperatures. This is largely because the total photoionization rate is approximately the same in all of our atmosphere models, meaning that the amount of energy in the non-thermal electron spectrum is also approximately constant. The volumetric heating rate for electrons is lower in simulations with higher CO₂ abundances because the electron densities are lower, but the heating rate per electron is in fact higher.

The lower panel in Fig. 12 shows the mixing ratios of several important species at the exobase as a function of CO₂ abundance, including the total ion mixing ratio. The exobase composition changes significantly, especially for N₂ and O₂; although O remains the most abundant species at the exobase in all simulations, in the 1000× simulation, the N₂ mixing ratio is similar to that of O. The changes in the exobase N₂ and O₂ mixing ratios are primarily because of the change in the exobase altitude. As the atmosphere cools and the exobase moves to lower altitudes, the molecular diffusion rates and the distance between the homopause and the exobase are reduced, giving molecular diffusion less of a chance to separate the heavier and lighter species.

Fig. 11 Neutral temperature structures for our simulations with different base mixing ratios of CO2 (upper-panel) and the total atmospheric heating rates due to the various heating processes as a function of CO2 mixing ratio (lower-panel). In the upper panel, the lines (which stop at the exobase) are simulations where the base CO2 density is 0.01, 0.1, 1, 10, 100, and 1000 times the value for the current Earth, which has a mixing ratio at 65 km of 4 × 10−4.

Fig. 12 Temperatures (upper panel), altitudes (middle panel), and mixing ratios of several important species of the exobase (lower panel) as a function of the CO2 mixing ratio at 65 km. In the lower panel, the black line gives the total ion mixing ratio (i.e. the ionization fraction).

4.2 The evolution of Earth’s upper atmosphere

In this section, we explore the responses of the upper atmospheres of Earth to the evolving XUV spectrum of the Sun between 3 Gyr in the past and 2.5 Gyr in the future. We do not consider the effects of the evolving lower atmospheric composition, which we will study in future work. We use the XUV spectra from Claire et al. (2012), who produced solar spectra as a function of age for all wavelengths that are of interest to us. We note that at young ages, the XUV spectra of solar mass stars are not unique functions of age given that they follow different rotational and activity evolution tracks (Johnstone et al. 2015a), which can be very important for the evolution of a planet’s atmosphere (Johnstone et al. 2015b). Since we do not know how rapidly the Sun was rotating at young ages, we do not know its early XUV output (Tu et al. 2015). However, at the ages that we consider, the rotation rates of young solar mass stars have converged to unique age dependent values, and so the unique XUV evolution presented by Claire et al. (2012) is valid.

In Fig. 13, we show our results for the Earth. In the upper-left panel, the temperature structures for our models at several ages are shown; the lines in this panel can be seen as an evolutionary sequence from right to left, with the thermospheres becoming cooler and less extended as the Sun’s XUV spectrum decays. In all cases, the neutral and ion temperatures are approximately the same at all altitudes and the electron temperatures are higher. In the upper-right panel, the evolutions of the exobase temperatures are shown. Interestingly, the electron temperatures at the exobase are higher by approximately the same amount (i.e. ~500 K) at all ages. At 3 Gyr ago, our models suggest that the exobase neutral and electron temperatures were approximately 4600 and 5100 K respectively. Our models also suggest that going into the future, we should not expect a large change in the upper atmosphere of the Earth simply due to the Sun’s activity decay. This is consistent with the power-law dependence of the solar activity on age, which means that the fastest changes take place at young ages; in our models, the atmosphere changes the most between 3 and 2 Gyr ago, and all changes after that are relatively slow in comparison.

Our results for the temperature structures are broadly consistent with those calculated by Tian et al. (2008a). Since Claire et al. (2012) used a more realistic method for estimating the spectra of the Sun at higher activity levels, we should not expect any exact agreement between our models and those of Tian et al. (2008a) even for the same total input XUV flux. Our model for 3 Gyr in the past corresponds approximately to their model with input fluxes of 4.9 times the current solar value, and we get similar results. Although not studied in this paper, in models withan even more active Sun, we can also see the effects of adiabatic cooling that were found by Tian et al. (2008a).

The middle row of Fig. 13 shows the exobase altitude and composition. The largest change in the altitude of the exobase takes place between 3 and 2 Gyr ago, dropping from 2000 km to approximately 600 km in that time. This is mostly due to the decrease in the thermospheric temperature, but is also partly due to the change in the chemical composition of the gas. At the youngest age considered, we find that the chemical composition at the exobase contains much more N and much less N₂ than we find for the current Earth, though O is in all cases the dominant species. Interestingly, we do not find any change in the exobase mixing ratios of H and He, which stay at values of ~ 5 × 10⁻⁵ at all ages. The ionization fraction at the exobase also decreases with age by about an order of magnitude between 3 Gyr ago and now.

Since the atmospheres are never fully hydrodynamic in our simulations, the dominant mass loss rates are likely to be non-thermal processes. Calculating these requires the application of additional models, which we do not attempt in this paper. We can calculate from our models the Jeans escape rates, which are shown for several species in the lower panels of Fig. 13. For the current Earth, it is well known that the only species undergoing significant Jeans escape is H, due to its low molecular mass. As we go to the past however, the Jeans escape rates of He, O, and N increase significantly, and in our earliest model are within approximately two orders of magnitude of the escape rate of H. If we were to go further into the past, the Jeans escape rates of these species would become comparable to that of H; it is around this time that we would find the atmosphere becoming hydrodynamic and the effects of adiabatic cooling on the temperature structures becoming important.

It is interesting to consider how the various heating mechanisms evolve with the evolving solar spectrum. In Fig. 14, we show the evolutions of direct XUV heating, chemical heating, and photoelectron heating for three altitude ranges. The quantities are the volumetric heating rates integrated over all the relevant altitudes and normalised to the modern Earth values. For altitudes above 100 km, the total heating rates for all these processes decrease with age due to the decreasing solar magnetic activity. However, at lower altitudes, the heating rates in fact increase with age. This is because this heating comes almost entirely from O₃ absorption at wavelengths longer than 2000 Å, which is generated in the solar photosphere. Since the Sun’s bolometric luminosity is increasing as it ages, the heating in the lower atmosphere of the Earth also increases. We note that the heating in the lower regions of our model should in fact be lower at young ages due to the absence of O₂ in the atmosphere prior to the Great Oxidation Event. The largest change is seen for heating of thermal electrons by non-thermalphotoelectrons, which we find was around 14 times larger 3 Gyr ago than the modern value. This is a result of the higher photoionization rates leading to there being more energy in the photoelectron spectrum and the greater ionization fractions of the gas leading to a larger fraction of that energy being transferred to the thermal electrons through elastic collisions, as opposed to being transferred to the neutrals through inelastic collisions.

Fig. 13 Our simulation results for the response of the Earth’s atmosphere to the evolving XUV spectrum of the Sun. In the upper left panel, we show the temperature structures for several of these models, where the different colours are for different ages and the solid, dashed, and dotted lines show the neutral, ion, and electron temperatures respectively. In all cases, the lines end at the exobase. In the upper right, middle left, and middle right panels, we show the exobase temperatures, altitudes, and chemical compositions as functions of age. In the lower left and lower right panels, we show the evolution of Jeans escape for H, O, N, and He, with the difference between the two plots being the range on the y-axis. In each figure, the small circles show the exact locations of each simulation.

Fig. 14 Evolution of direct XUV heating (upper panel), chemical heating (middle panel), and photoelectron heating (lower panel) within three altitude ranges in the Earth’s atmosphere. The quantities plotted are the total heating rates integrated over the relevant volumes and normalised to the values for the modern Earth. For photoelectron heating, only the values for altitudesbetween 200 km and the exobase are shown since at lower altitudes this process is negligible.

5 Discussion

In this paper, we develop and validate a physical model for the thermal and chemical structures of planetary upper atmospheres. We use this model to estimate how the Earth’s upper atmosphere would look if the CO₂ abundance was varied by a large amount, and to explore the response of the upper atmosphere to the evolving XUV spectrum of the Sun. Increasing the CO₂ abundance causes the Earth’s upper atmosphere to be much cooler and contracted, presumably causing a dramatic decrease in atmospheric losses. This could also be important for the evolution of the Earth’s atmosphere since without some form of protection, such as the enhanced CO₂ abundances expected during earlier epochs (Zahnle et al. 2010), interactions with the solar wind could have stripped away the early Earth’s atmosphere (Lichtenegger et al. 2010). This is also important because, based on the examples of Venus, Mars, and likely the early Earth, CO₂ dominated atmospheres are likely common among terrestrial exoplanets.

Our physical model should be applicable to a range of situations with arbitrary atmospheric compositions and stellar input XUV and IR spectra. For this purpose, we have attempted to develop the model in such a way that, as much as possible, it is based on first-principles physics. This means that we have avoided adding free parameters in our model, and we rely as little as possible on parametrised scaling laws that have been developed for the solar system terrestrial planets. For example, our model does not rely on parametarised heating efficiencies, but instead it calculates the heating rates for each of the contributing physical mechanisms explicitly. However, our model still contains several weaknesses; for example, our treatment IR cooling should be improved with more detailed treatment of radiation transfer and the detailed interactions between the IR field and the various coolants.

The main reason that we are interested in the thermal and chemical structures of the atmospheres of terrestrial planets is their importance for atmospheric losses into space. Even just for the case of the Earth, the atmosphere has radically changed many times, especially in the first billion years after the Earth’s formation. Models such as ours are essential tools for understanding this diverse set of atmospheric conditions. Estimating atmospheric loss rates atmospheres would require additional detailed modelling of the planetary exospheres (e.g. Kislyakova et al. 2013) and of the interactions of the ionosphere with the planet’s magnetic field (e.g. Glocer et al. 2009), which we do not attempt in this paper. In future work, we will combine the model developed in this paper with other such models to gain a more complete understanding of how losses to space influence atmospheric evolution.

Acknowledgments

The authors thank the anonymous referee for useful comments on our manuscript. This study was carried out with thesupport by the FWF NFN project S11601-N16 “Pathways to Habitability: From Disk to Active Stars, Planets and Life” and the related subprojects S11604-N16, and S11607-N16. We also acknowledge the support by grant 16-52-14006 of the Russian Fund of Basic Research, and the Austrian Science Foundation (FWF) projects I2939-N27, S11606-N16.

Appendix A The Kompot code

To solve the physical model described above, we have developed The Kompot Code, which calculates the 1D thermal and chemical structure of a planet’s atmosphere. The underlying equations for the physical model are described in Sect. 2, and the numerical methods used to solve these equations are described in the following appendices. The code is designed to be very flexible, such that it can be applied to a wide range of atmospheres, and the underlying physics can be easily modified and improved. The physical input parameters into the code are the temperature and species densities at the lower boundary, and the stellar XUV and IR spectra at the exobase. It is written in Fortran and Python and has been parallelised using OpenMP. The Kompot Code will be made publicly available in the near future and will be obtainable by contacting the authors.

Evolving the state of the atmosphere forward in time using the system of equations described in Sect. 2 is not trivial since they contain many terms that are different in form. For example, consider Eq. (1) in the following form:- $$\frac{\partial n_j}{\partial t}=-\frac{1}{r^2}\frac{\partial \left(r^2{n}_{j}v\right)}{\partial r}-\frac{1}{r^2}\frac{\partial \left(r^2{\Phi }_{\text{d},j}\right)}{\partial r}+S_j.$$(A.1)

The first term on the RHS is for advection and contains a ∂n_j ∕∂r term. The second term is for diffusion and contains multiple ∂²n_j∕∂r² terms. The third term describes the chemical sources and is a stiff system of many ordinary differential equations. No one numerical method is ideally suited to solve all of these problems simultaneously. We instead use operator splitting to solve each one sequentially using different numerical methods.

The simulation takes place on a finite static spatial grid of cells and we calculate all quantities at cell centres. For accuracy and efficiency we place our cells in such a way that the cell spacing increases linearly with altitude. The code only considers the cells below the exobase, the location of which varies within the simulation.

The initial conditions of a simulation should be irrelevant for the final result, but this is only the case when they are not too unrealistic. We have three ways to calculate the initial conditions. Firstly, we can start the simulation assuming hydrostatic equilibrium, a uniform temperature, and uniform mixing ratios for all species. For numerical reasons, it is often necessary to evolve the initial chemical structure of the atmosphere by an arbitrary amount of time before starting the simulation, mostly to avoid starting with no ions and electrons in the gas. Secondly, we can start many simulations simply using the atmospheric structure from the results of previous simulations. Thirdly, in fully hydrodynamic simulations, we start the simulation assuming an isothermal Parker wind structure and uniform chemical composition.

We evolve the atmosphere forward in time by performing a large number of timesteps, with length Δ t. We determine Δt by multiplying the minimum time taken for a sound wave to cross one grid cell by a fixed number, called the Courant number, C. Unless we are applying the full hydrodynamic method given in Appendix B, our simulations are not strictly constrained by the Courant Friedrichs Lewy condition that C should be less than unity, though we often assume C ~ 1 to ensure numerical stability. At the beginning of a timestep, we first update the XUV, IR, and non-thermal electron spectra at each grid cell; these are then assumed to be constant for the entire timestep. We then update the atmospheric properties by applying each of the physical mechanisms successively in the following order: hydrodynamics, chemistry, diffusion, heating/cooling, energy exchange, conduction. At the end of each timestep, we then recalculate the exobase location.

In orderto make the simulations computationally more efficient, we generally do not do every part of the above set of steps in every timestep. It is not necessary to recalculate the XUV and photoelectron spectra and the IR flux every timestep. Instead, we update these quantities every 100 timesteps. Similarly, we find that applying the chemical network to evolve the species densities is very time consuming and we get almost identical results by only making this step every 100 timesteps. When we do evolve the chemistry, we do so for the entire 100 timesteps since it was last evolved. In most applications of our model, we are interested in the steady-state atmospheric conditions that result in constant input parameters. To calculate this, we perform successive timesteps starting from the initial conditions until the atmosphere reaches a steady state. Our model could also be used to simulate the time evolution of the atmosphere in response to changing input conditions (e.g. a stellar flare).

Appendix B Hydrodynamics solver

In this appendix, we give the algorithm used for solving the full set of time-dependent hydrodynamic equations. The algorithm used is an explicit solver that is 2nd order accurate in time and 3rd order accurate in space. Although not used in this paper, we describe the solver here because it will be used in future studies with this model.

The hydrodynamic equations in spherical coordinates can be written $$\begin{array}{lll}\hfill & \hfill & \frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v\right)}{\partial r}=-\frac{2\rho v}{r},\hfill \\ \hfill & \hfill & \frac{\partial (\rho v)}{\partial t}+\frac{\partial \left(\rho vv\right)}{\partial r}=-\frac{\partial p}{\partial r}-\frac{2\rho v^2}{r}-\rho g,\hfill \\ \hfill & \hfill & \frac{\partial e_{\text{n}}}{\partial t}+\frac{\partial \left(e_{\text{n}}v\right)}{\partial r}=-\frac{\partial \left(p_{\text{n}}v\right)}{\partial r}-\frac{2v\left(e_{\text{n}}+p_{\text{n}}\right)}{r}-{\rho }_{\text{n}}vg,\hfill \\ \hfill & \hfill & \frac{\partial e_{\text{i}}}{\partial t}+\frac{\partial \left(e_{\text{i}}v\right)}{\partial r}=-\frac{\partial \left(p_{\text{i}}v\right)}{\partial r}-\frac{2v\left(e_{\text{i}}+p_{\text{i}}\right)}{r}-{\rho }_{\text{i}}vg,\hfill \\ \hfill & \hfill & \frac{\partial e_{\text{e}}}{\partial t}+\frac{\partial \left(e_{\text{e}}v\right)}{\partial r}=-\frac{\partial \left(p_{\text{e}}v\right)}{\partial r}-\frac{2v\left(e_{\text{e}}+p_{\text{e}}\right)}{r}-{\rho }_{\text{e}}vg,\hfill \end{array}$$

The LHSs have the forms of pure advection equations and the RHSs contain source terms for pressure (terms involving ∂∕∂r), the spherical geometry (terms involving 1∕r), and gravity (terms involving g). Writing the equations in this way simplifies the problem significantly since we can then solve the advection and sources separately.

The above equations in vector form are $$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial r}=S,$$(B.6)

where $$U=\left(\begin{array}{c}\rho \\ \rho v\\ e_{\text{n}}\\ e_{\text{i}}\\ e_{\text{e}}\end{array}\right),F=vU,S=\left(\begin{array}{c}-\frac{2\rho v}{r}\\ -\frac{\partial p}{\partial r}-\frac{2\rho v^2}{r}-\rho g\\ -\frac{\partial \left(p_{\text{n}}v\right)}{\partial r}-\frac{2v\left(e_{\text{n}}+p_{\text{n}}\right)}{r}-{\rho }_{\text{n}}vg\\ -\frac{\partial \left(p_{\text{i}}v\right)}{\partial r}-\frac{2v\left(e_{\text{i}}+p_{\text{i}}\right)}{r}-{\rho }_{\text{i}}vg\\ -\frac{\partial \left(p_{\text{e}}v\right)}{\partial r}-\frac{2v\left(e_{\text{e}}+p_{\text{e}}\right)}{r}-{\rho }_{\text{e}}vg\end{array}\right).$$(B.7)

Here, U are the conserved quantities that we want to update, F are the advection fluxes of these quantities, and S are the sources. In the description below, we use the subscript j to refer to the cell index and the superscript n to refer to the timestep number, such that $U_j^n$ is the value of U at the jth cell and the nth timestep. The subscript j + 1∕2 is used to refer to quantities at the boundary between the jth and (j + 1)th cells at radius r_j+1∕2 = (r_j+1 − r_j)∕2.

For the time integration, we split these equations into advection and source steps using Strang splitting. This means that we first evolve U by half a timestep due to the sources, then evolve the updated values of U a full timestep due to advection, and finally evolve U again a halftimestep due to the sources⁵. For each one of these updates, we use the 2nd order Total Variation Diminishing (TVD) Runge–Kutte scheme given by Gottlieb & Shu (1998), given by $$\begin{array}{lll}\hfill & \hfill & U^{(1)}=U^n+\Delta tf(U^n),\hfill \\ \hfill & \hfill & U^{n+1}=\frac{1}{2}\left(U^n+U^{(1)}\right)+\frac{1}{2}\Delta tf(U^{(1)}),\hfill \end{array}$$

where f(Uⁿ) is ∂U ∕∂t calculated from Uⁿ. The first step is the standard Forward Euler method, and the second step improves the approximation. When doing the two source term updates, Δ t above should be replaced with Δt∕2.

For the advection part (i.e. S = 0), Eq. (B.6) can be expressed in terms of the boundary fluxes as $$f(U_j^n)=-{\left(\frac{\partial F}{\partial r}\right)}_j^n=-\frac{F_{j+\frac{1}{2}}^n-F_{j-\frac{1}{2}}^n}{\Delta r_j},$$(B.10)

where Δr_j = r_j+1∕2 − r_j−1∕2 is the cell width. To calculate the cell boundary fluxes, we use the MUSCL approach using the high-resolution TVD Lax–Friedrichs numerical flux and the minmod slope limiter (see Sects. 3.5.1 and 3.5.4 of Yee 1989). In the following discussion, we do not write the superscript n and all quantities should be assumed to be from the nth timestep. The basic idea of the MUSCL approach is instead of calculating the flux $F_{j+\frac{1}{2}}$ from the adjacent cell centre values of U (i.e. U_j and U_j+1), we calculate it from left and right values of $U_{j+\frac{1}{2}}$, given by $$\begin{array}{lll}\hfill & \hfill & U_{j+\frac{1}{2}}^{\text{L}}=U_j+\frac{1}{4}\left[\left(1-\eta \right){\tilde{\tilde{\Delta }}}_{j-\frac{1}{2}}+\left(1+\eta \right){\tilde{\Delta }}_{j+\frac{1}{2}}\right],\hfill \\ \hfill & \hfill & U_{j+\frac{1}{2}}^{\text{R}}=U_{j+1}-\frac{1}{4}\left[\left(1-\eta \right){\tilde{\Delta }}_{j+\frac{3}{2}}+\left(1+\eta \right){\tilde{\tilde{\Delta }}}_{j+\frac{1}{2}}\right],\hfill \end{array}$$

where $$\begin{array}{lll}\hfill & \hfill & {\tilde{\Delta }}_{j+\frac{1}{2}}=\text{minmod}({\Delta }_{j+\frac{1}{2}},\omega {\Delta }_{j-\frac{1}{2}}),\hfill \\ \hfill & \hfill & {\tilde{\tilde{\Delta }}}_{j+\frac{1}{2}}=\text{minmod}({\Delta }_{j+\frac{1}{2}},\omega {\Delta }_{j+\frac{3}{2}}),\hfill \end{array}$$

and Δ_j+1∕2 = U_j+1 −U_j. The quantities η and ω are discussed below. The minmod slope limiter is given by $$\text{minmod}(x,\omega y)=\frac{\left|x\right|}{x}\max \left[0,\min \left(\left|x\right|,\frac{\omega y\left|x\right|}{x}\right)\right].$$(B.15)

In the MUSCL scheme, the TVD Lax–Friedrichs numerical flux is $$\begin{array}{lll}{F}_{j+\frac{1}{2}}\hfill & =\hfill & \frac{1}{2}\left[F\left(U_{j+\frac{1}{2}}^{\text{L}}\right)+F\left(U_{j+\frac{1}{2}}^{\text{R}}\right)\right]\hfill \\ \hfill & \hfill & -\frac{\Delta t}{2\Delta r_{j+\frac{1}{2}}}\left(U_{j+\frac{1}{2}}^{\text{R}}-U_{j+\frac{1}{2}}^{\text{L}}\right),\hfill \end{array}$$(B.16)

where $\Delta r_{j+\frac{1}{2}}=r_{j+1}-r_j$. To calculate the fluxes on the RHS of this equation, we use $$\begin{array}{lll}\hfill & \hfill & F\left(U_{j+\frac{1}{2}}^{\text{L}}\right)=\frac{{\left(\rho v\right)}_{j+\frac{1}{2}}^{\text{L}}}{{\rho }_{j+\frac{1}{2}}^{\text{L}}}{U}_{j+\frac{1}{2}}^{\text{L}},\hfill \\ \hfill & \hfill & F\left(U_{j+\frac{1}{2}}^{\text{R}}\right)=\frac{{\left(\rho v\right)}_{j+\frac{1}{2}}^{\text{R}}}{{\rho }_{j+\frac{1}{2}}^{\text{R}}}{U}_{j+\frac{1}{2}}^{\text{R}}.\hfill \end{array}$$

Here, ${\left(\rho v\right)}_{j+\frac{1}{2}}^{\text{L}}$ and ${\rho }_{j+\frac{1}{2}}^{\text{L}}$ are elements of $U_{j+\frac{1}{2}}^{\text{L}}$ and ${\left(\rho v\right)}_{j+\frac{1}{2}}^{\text{R}}$ and ${\rho }_{j+\frac{1}{2}}^{\text{R}}$ are elements of $U_{j+\frac{1}{2}}^{\text{R}}$. The value of η chosen determines the spatial order of accuracy of the scheme; we take η = 1∕3, making the scheme third-order accurate in space with an upwind bias. Given that value of η, the value of ω is an adjustable parameter with a maximum of 4; we assume ω = 3.5.

The equation for the source term updates is $$f(U_j^n)=S_j^n,$$(B.19)

where S is given by Eq. (B.7). The only assumption that needs to be made to discreetise S is for the pressure source terms since they contain spatial derivatives. For these terms, we simply use central differencing.

We have performed several tests of our implementation of this algorithm, including using it to calculate the structure of a 1D isothermal stellar wind with a known analytic solution (Parker 1958). In this test, our calculations match the analytic solutions exactly. We have tested our hydrodynamics solver using the Versatile Advection Code (VAC; Tóth 1996) and find that it performs well in the standard Sod shock tube test, though VAC is less diffusive at low spatial resolution. Finally, we are able to reproduce well the hydrodynamic atmosphere simulations of Johnstone et al. (2015b), also performed using VAC.

We must also evolve the densities of individual species due to advection. To do this, we convert the cell boundary fluxes for mass into cell boundary fluxes for individual species using a simple upwind approximation for the mixing ratios of individual species at the cell boundaries. Specifically, the mixing ratios for each species are assumed to be equal to those at the centre for the cell that the mass is flowing from, meaning that if the flux is positive at the j + 1∕2 boundary, then the mixing ratio at the jth cell is taken.For the kth species, the flux is $$F_{k,j+1/2}=\{\begin{array}{ll}{n}_{k,j}{F}_{\rho ,j+1/2}/{\rho }_{k,j+1/2},\hfill & \text{if}{F}_{\rho ,j+1/2}>0,\hfill \\ 0,\hfill & \text{if}{F}_{\rho ,j+1/2}=0,\hfill \\ n_{k,j+1}{F}_{\rho ,j+1/2}/{\rho }_{k,j+1/2},\hfill & \text{if}{F}_{\rho ,j+1/2}<0.\hfill \end{array}$$(B.20)

With these fluxes, we then calculate the update for the densities in each cell using a simple forward Euler method for the time integration, such that $$n_{k,j}^{n+1}=n_{k,j}^n-\frac{\Delta t}{\Delta r_j}\left(F_{k,j+\frac{1}{2}}-F_{k,j-\frac{1}{2}}\right).$$(B.21)

After this update, small inconsistencies between the values of ρ and the values of n_k at each cell are corrected for by scaling the n_k values such that ∑_km_kn_k = ρ.

Appendix C Tridiagonal matrix algorithm

The tridiagonal matrix algorithm is a commonly used algorithm and detailed descriptions of the method can be found in many textbooks. Since we apply this solver for several different parts of the model, it is necessary to repeat the main steps in the algorithm here for the explanations in the following appendices to make sense. The form of the algorithm presented here is adapted from Bodenheimer et al. (2007).

The aim of the algorithm is to solve the tridiagonal systems of equations, which can be written as a system of simultaneous equations, each given by $$a_j{x}_{j-1}+b_j{x}_j+c_j{x}_{j+1}=-d_j.$$(C.1)

In our solvers, we have one of these equations for each grid cell. In the tridiagonal matrix algorithm, the coefficients a_j, b_j, c_j, and d_j are known, and the aim is to calculate x_j, which represents the physical quantity of interest. The aim in the next three appendices is to derive expressions for these coefficients, and H_j and Y_j discussed below, for the different physical mechanisms.

To solvethis system of equations, consider the equation $$x_j=H_{j-1}{x}_{j-1}+Y_{j-1}.$$(C.2)

Firstly, the values of H_j and Y_j need to be calculated for each value of j except for j = J. Assuming we know H_J−1 and Y_J−1 which are derived separately for each problem, the other values are calculated by iterating downwards through the grid using $$H_{j-1}=-\frac{a_j}{b_j+c_j{H}_j}$$(C.3)

and $$Y_{j-1}=-\frac{d_j+c_i{Y}_i}{b_j+c_j{H}_j}.$$(C.4)

With H_j and Y_j known, the values of x_j is calculated at each cell by iterating upwards using Eq. (C.2) and assuming that the value of x₁ is known in advance (this is always the case since we use fixed lower boundary values for all quantities).

Appendix D Solver for semi-static hydrodynamic equations

We discuss in this appendix first the method used to solve the energy equations and then the method used to solve the equations for ρ and v when solving the hydrodynamics using the semi-static method presented in Sect. 2.2. It is convenient to update the energy each timestep using an implicit method since this avoids the restrictively short timestep sizes needed in the explicit schemes. The hydrodynamic part of the energy equation including gravity is $$\frac{\partial e}{\partial t}=-\frac{1}{r^2}\frac{\partial \left[r^2{\Phi }_e\right]}{\partial r}-\rho vg,$$(D.1)

where ${\Phi }_e=v\left(e+p\right)$. For this section, we do not write the subscripts n, i, and e for the neutral, ion, and electron components; the method presented here is used to solve the energy equations for each component separately. Using the Crank–Nicolson method for time discreetization, and the cell boundary values to discreetise cell centered spatial derivatives, we get $$\begin{array}{lll}\frac{e_j^{n+2}-e_j^n}{\Delta t}\hfill & =\hfill & -\frac{1}{2r_j^2}\frac{\left(r_{j+\frac{1}{2}}^2{\Phi }_{e,j+\frac{1}{2}}^{n+1}-r_{j-\frac{1}{2}}^2{\Phi }_{e,j-\frac{1}{2}}^{n+1}\right)}{\Delta r_j}\hfill \\ \hfill & \hfill & -\frac{1}{2r_j^2}\frac{\left(r_{j+\frac{1}{2}}^2{\Phi }_{e,j+\frac{1}{2}}^n-r_{j-\frac{1}{2}}^2{\Phi }_{e,j-\frac{1}{2}}^n\right)}{\Delta r_j}-{\rho }_j{v}_j{g}_j.\hfill \end{array}$$(D.2)

This can be rewritten $$\begin{array}{lll}{e}_j^{n+2}+k_j\left(r_{j+\frac{1}{2}}^2{\Phi }_{e,j+\frac{1}{2}}^{n+1}-r_{j-\frac{1}{2}}^2{\Phi }_{e,j-\frac{1}{2}}^{n+1}\right)\hfill & =\hfill & e_j^n-k_j\hfill \\ \hfill & \hfill & \times \left(r_{j+\frac{1}{2}}^2{\Phi }_{e,j+\frac{1}{2}}^n-r_{j-\frac{1}{2}}^2{\Phi }_{e,j-\frac{1}{2}}^n\right)\hfill \\ \hfill & \hfill & -{\rho }_j{v}_j{g}_j\Delta t,\hfill \end{array}$$(D.3)

where $$k_j=\frac{\Delta t}{2r_j^2\Delta r_j}.$$(D.4)

To get Φ_e at cell boundaries, we assume that the cell boundary values of spatially variable quantities are the average of the cell centre values; for example $e_{j+\frac{1}{2}}=(e_{j+1}+e_j)/2$. Therefore, the cell boundary fluxes is given by $$\begin{array}{lll}{\Phi }_{e,j+\frac{1}{2}}^n\hfill & =\hfill & \frac{1}{2}{\gamma }_{j+\frac{1}{2}}{v}_{j+\frac{1}{2}}{e}_j^n+\frac{1}{2}{\gamma }_{j+\frac{1}{2}}{v}_{j+\frac{1}{2}}{e}_{j+1}^n\hfill \\ \hfill & \hfill & -\frac{1}{2}({\gamma }_{j+\frac{1}{2}}-1){\rho }_{j+\frac{1}{2}}{v}_{j+\frac{1}{2}}^3.\hfill \end{array}$$(D.5)

where we have used ${\Phi }_e=v(e+p)=\gamma ev-\frac{1}{2}(\gamma -1)\rho v^3$. The equations for ${\Phi }_{e,j-\frac{1}{2}}^n$, ${\Phi }_{e,j+\frac{1}{2}}^{n+1}$, and ${\Phi }_{e,j-\frac{1}{2}}^{n+1}$ have the same form and can be obtained with the appropriate substitutions of the subscripts and superscripts.