Principles of atmospheric modeling with Exo_k

Author:

Jérémy Leconte (LAB/CNRS/Univ. Bordeaux)

Exo_k includes a class dedicated to model self consistently the evolution of planetary atmospheres: the Atm_evolution class. The heart of this model is the native Exo_k radiative transfer core so that any available radiative data type and format can be used in the evolution module. In addition to radiation, the atmospheric model includes the effect of dry convection, turbulent diffusion, and moist processes (condensation, moist convection, and precipitations) for any number of species.

As in a number of problems of interest, only the equilibrium state of the atmosphere is of interest, a particular attention has been devoted to the development of schemes to accelerate the convergence toward such steady-state.

This section focuses on the physical description of the model. To learn more about the way you can use the library to model atmospheres, have a look at the tutorial.

Basic structure of the atmosphere: the Atm class.

In Exo_k, the state of the atmosphere at any given moment is described by the Atm class that handles, among other things, the radiative transfer calculations. As can be seen on the schematic, the atmosphere is composed of \(N_\mathrm{lay}\) layers at temperature \(T_{\mathrm{n}}\) where \(\mathrm{n}\) goes from 0 at the top of atmosphere to \(N_\mathrm{lay}-1\) for the layer just above the surface. The layer \(\mathrm{n}\) is bounded at the top by a level interface at pressure \(p_\mathrm{lev,\mathrm{n}}\) and at the bottom by a level at \(p_\mathrm{lev,\mathrm{n}+1}\). The pressure at the middle of the layer is called \(p_\mathrm{lay,n}\). For sake of conciseness, when possible, we will directly refer to the temperature and pressure vectors defined as

\[ \begin{align}\begin{aligned}\begin{aligned}\\\hat{T}= \{T_{\mathrm{n}}\}_{\mathrm{n}\in\{0,\,N_\mathrm{lay}-1\}}\end{aligned}\end{aligned}\end{align} \]

\[\begin{aligned} \hat{P}_\mathrm{lay}= \{p_\mathrm{lay}\}_{\mathrm{n}\in\{0,\,N_\mathrm{lay}-1\}}.\end{aligned}\]

For quantities defined at levels, the vectors have one additional element, such as

\[\begin{aligned} &\hat{P}_\mathrm{lev}= \{p_\mathrm{lev}\}_{\mathrm{n}\in\{0,\,N_\mathrm{lay}\}}.\end{aligned}\]

The same convention apply for other vectors.

Schematic of the organisation of model layers and levels in a case with \(N_\mathrm{lay}=4\).

The mass per unit area inside each layer is given by

\[\begin{aligned} \hat{M}= \frac{\Delta\hat{P}_\mathrm{lev}}{g}, \label{masses}\end{aligned}\]

where \(g\) is the gravity of the planet [1] and the \(\Delta\) operator is defined such that the \(\mathrm{n}\)-th component of \(\Delta\hat{X}\) is equal to \(X_{\mathrm{n}+1}-X_\mathrm{n}\). \(\Delta\hat{X}\) thus has one less element. We highlight the fact that, in our convention, layers are counted from the top down. The \(\Delta\) operator thus yields the value of a quantity in a given layer minus the quantity in the layer above it. As a result, the layer masses in Eq. ([masses]) are positive.

Radiative transfer

The radiative transfer is computed using a 2-stream approximation of Toon et al. (1989) using \(N_\mathrm{lay}-1\) radiative layers. The \(\mathrm{n}\)-th radiative layer goes from the middle of the \(\mathrm{n}\)-th layer (\(p_\mathrm{lay,n}\)) to the middle of the layer below it (\(p_\mathrm{lay,\mathrm{n}+1}\)). This choice is motivated by the fact that we need to know the source function (i.e. the temperature) at the interfaces between the radiative layers for the algorithm to be stable.

The opacities of the radiative layers are computed using the Exo_k opacity library and can use all the types of opacity sources available through it (correlated-\(k\), cross sections, collision-induced absorptions, Rayleigh scattering, aerosol Mie scattering, etc.). They are computed at level surfaces with temperatures and volumic concentrations of the various species composing the atmosphere (\(x\)) computed as shown in Fig. [fig:levels].

Because we have decided to take \(p_\mathrm{lay,N_\mathrm{lay}-1} \equiv p_\mathrm{lev,N_\mathrm{lay}}\), the bottom of the last radiative layer is effectively the surface and we treat it as an opaque surface with an albedo, \(A_\mathrm{s}\), that depends on the wavenumber (with an emissivity equal to \(\epsilon_\mathrm{s}=1-A_\mathrm{s}\)). At the top, we set the diffuse downwelling flux to be equal to the incoming flux set by the user (e.g. the average incoming stellar flux impinging on the planet). This means that all sources of opacity between the top of the model and the \(p=0\) level are disregarded. This can lead to an overestimated cooling to space of the top layer.

The net radiative fluxes, \(\hat{F}_\mathrm{net}\), are then computed at the level interfaces (i.e. at the middle of the radiative layers) using Toon et al. (1989). Finally, the heating rates (in power per unit mass of air, or W/kg) are given by

\[\begin{aligned} \hat{H}^\mathrm{rad}=\Delta\hat{F}_\mathrm{net}/ \hat{M}.\end{aligned}\]

The internal flux coming form the planet’s interior, if any, is added to the budget of the bottom layer.

Time integration

Once we are able to compute heating rates, we can compute the thermal evolution of the atmosphere. This time intégration is handled by the Atm_evolution class.

The conservation of energy states that the time derivative of the temperature (\(\dot{\hat{T}}\)) is given by

\[\begin{aligned} c_p\dot{ \hat{T}} = \hat{H},\end{aligned}\]

where \(\hat{H}\) are the heating rates expressed in power per unit mass of air. In the model, because \(c_p\) is considered constant throughout the atmosphere, it is convenient to define a reduced time following \(t^{*}\equiv t/c_p\) so that the equation reads

\[\begin{aligned} \hat{T}^\prime = \hat{H}, \label{time_derivative}\end{aligned}\]

where \(^\prime\) will always refer to a derivative with respect to this new, reduced time. In the following, any time duration or timescale whose symbol is accompanied by a \(^*\) will be implicitly understood as using this reduced time unit as well.

For each physical process that is implemented in the model, the relevant parametrization can be seen as an operator that provides a heating rate in each layer considering the current thermal state of the atmosphere, denoted by \(\hat{T}^{t^{*}}\). The thermal state after a time \(\Delta t^{*}\) is thus

\[\begin{aligned} \hat{T}^{t^{*}+\Delta t^{*}}= \hat{T}^{t^{*}}+ \Delta t^{*}\cdot \hat{H}(\hat{T}^{t^{*}}) . \label{T_operator}\end{aligned}\]

For sake of compactness, for each physical process, we can define a new, related operator so that

\[\begin{aligned} \Xi^{\Delta t^{*}}(\hat{T}^{t^{*}})\equiv \hat{T}^{t^{*}}+ \Delta t^{*}\cdot \hat{H}(\hat{T}^{t^{*}}) \Rightarrow \hat{T}^{t^{*}+\Delta t^{*}}= \Xi^{\Delta t^{*}}(\hat{T}^{t^{*}}) .\end{aligned}\]

At this stage, it is important to recognize that there are other variables of interest that are needed to compute the heating rates and that need to be integrated in time. Let us call \(\hat{\phi}\) the vector that describes the values of all the quantities of interest (including temperature and pressure) in all the layers of the model and generalize the \(\Xi^{\Delta t^{*}}\) operator so that

\[\begin{aligned} \hat{\phi}^{t^{*}+\Delta t^{*}}= \Xi^{\Delta t^{*}}(\hat{\phi}^{t^{*}}),\end{aligned}\]

where the thermal part will take a form similar to Eq. ([T_operator])

Instead of just adding the various contributions of the various physical parametrization evaluated on the initial state of a given timestep, we find that the algorithm is much more stable if the integration is done following

\[\begin{aligned} \hat{\phi}^{t^{*}+\Delta t^{*}}= \Xi^{\Delta t^{*}}_\mathrm{rain}\left(\Xi^{\Delta t^{*}}_\mathrm{cond}\left(\Xi^{\Delta t^{*}}_\mathrm{madj}\left(\Xi^{\Delta t^{*}}_\mathrm{conv}\left(\Xi^{\Delta t^{*}}_\mathrm{rad}(\hat{\phi}^{t^{*}})\right)\right)\right)\right),\end{aligned}\]

where \(\mathrm{rad}\), \(\mathrm{conv}\), \(\mathrm{madj}\), \(\mathrm{cond}\), and \(\mathrm{rain}\) stand for radiation, dry convection, moist convection, large-scale condensation, and rain. Each of these processes will be discussed in detailed later on. In other words, the atmospheric state at the end of each process is used as the initial state for the following physical process.

For future reference, let us clarify that, even though the heating rates due to the various processes have not been computed using the same thermal state, the thermal evolution still verifies

\[\begin{aligned} \hat{T}^{t^{*}+\Delta t^{*}}= \hat{T}^{t^{*}}+ \Delta t^{*}\left(\hat{H}^\mathrm{rad} + \hat{H}^\mathrm{conv} + \hat{H}^\mathrm{madj} + \hat{H}^\mathrm{cond} + \hat{H}^\mathrm{rain} \right).\end{aligned}\]

Physical processes

Here, we give a brief description of the various physical processes included in the model.

Dry convective adjustment

The dry convective adjustment scheme looks for convectively unstable regions in the atmosphere and brings them back to neutral stability. In general, these regions are identified by layers where the potential temperature, \(\hat{\theta}= \hat{T}(p_0/\hat{P}_\mathrm{lay})^{R/c_p}\), decreases upward. Here \(p_0\) is a constant reference pressure and \(R\) is the specific gas constant for the dry air.

However, this criterion is insufficient when there can be large variations of the mean molar mass, \(\hat{M}\), of the atmosphere that can affect the density of the gas. We thus use the virtual potential temperature, \(\hat{\theta}_\mathrm{v}= \hat{\theta}(M_\mathrm{a}/\hat{M})\), where \(M_\mathrm{a}\) is the mean molar mass of the dry air. Unstable layers are defined by \(\Delta\hat{\theta}_\mathrm{v}>0\).

The potential temperature and composition of unstable layers are fully mixed over a single timestep while conserving total mass and enthalpy.

Moist convective adjustment

Our scheme to handle moist convection is similar to the one used in Leconte et al. (2013) but has been generalized to any number of condensing species, although each species is treated separately.

To account for the possible convection inhibition due to molar mass effects, we follow the method from Leconte et al. (2017) and suppress the convective adjustment in any layer where the vapor mixing ratio for the condensing vapor exceeds the critical specific concentration given by

\[\begin{aligned} q_\mathrm{cri}\equiv \frac{RT}{\left(M_\mathrm{v}-M_\mathrm{a}\right) L}.\end{aligned}\]

Condensation

Condensation of vapor can also happen in a layer in absence of convective processes, when there is diabatic cooling, for example. To account for this, at each timestep, this parametrization brings any supersaturated layer back to vapor equilibrium. Subsaturated layers can also evaporate condensates if they are present. This saturation adjustment is performed iteratively until equilibrium conditions are found at constant moist enthalpy.

Rains

In the current version, all condensates precipitate instantaneously upon condensation, always leaving a cloud-free atmosphere. However, as we are interested in deep atmospheres where precipitations are unlikely to reach the surface (if it exists), we implemented a simple scheme to reevaporate falling precipitations. We start from the model top layer and collect all precipitations downward. Whenever an unsaturated layer is met, a fraction \(f\) of the condensate that would need to be evaporated to saturate the layer is effectively evaporated. Precipitation are fully evaporated when they reach a layer where the temperature is above the boiling temperature of the falling species. All remaining precipitations, if any, when the surface layer is reached are added to the surface.

Numerical acceleration

Computation of the adaptive timestep

As radiative transfer is usually the most expensive part of a 1D model, extra care has been taken to compute it as seldom as possible. As a first step in that direction, Exo_k has an adaptive timestep that is based on the radiative timescale of the atmosphere. This radiative timescale is computed as follows. We start by saying that, around a given thermal state of the atmosphere, \(\hat{T}_\mathbb{K}\), the heating rates can be linearized through

\[\begin{aligned} \hat{H}^\mathrm{rad}(\hat{T}_\mathbb{K}+\delta \hat{T}) = \hat{H}^\mathrm{rad}(\hat{T}_\mathbb{K}) + \mathbb{K}\cdot \delta \hat{T}, \label{jacobian}\end{aligned}\]

where \(\mathbb{K}\) is the Jacobian matrix of the heating rates (of dimension \(N_\mathrm{lay}\times N_\mathrm{lay}\)).

This matrix tells us how heating in a layer is related to a temperature change anywhere in our atmosphere. The biggest terms are the diagonal ones, which are usually negative. This shows that any layer whose temperature is increased will tend to emit more and cool in response. The terms directly above and below the diagonal – the coupling terms between adjacent layers – are usually positive but smaller in magnitude. Other off-diagonal terms, which are long range couplings, are much smaller, especially in optically thick atmospheres.

Now, let us consider the radiative evolution of the system around a thermal state in equilibrium — meaning that \(\hat{H}^\mathrm{rad}(\hat{T}_\mathrm{eq})=\vec{0}\). We can insert Eq. ([jacobian]) in Eq. ([time_derivative]), which yields

\[\begin{aligned} \delta \hat{T}^\prime = \mathbb{K}\cdot \delta \hat{T}.\end{aligned}\]

Keeping only the most important terms (the diagonal ones), the layers decouple and the solution reads

\[\begin{aligned} \delta \hat{T}= \delta \hat{T}_0 \cdot e^{- \mathrm{Diag} (\mathbb{K}) t^{*}}.\end{aligned}\]

So we see that in each layer, an initial perturbation will be radiated away on a reduced timescale equal to \(\tau^{*}=1/ \mathrm{Diag} (\mathbb{K})\) (The timescale in physical time is obtained with \(\tau=c_p\tau^{*}\)). In the baseline evolution, the smallest radiative timescale in our atmosphere is used as timestep to ensure that the radiative evolution is well sampled throughout the atmosphere. We find that this condition is sufficient to ensure a stable evolution in most cases, but the user can always specify smaller timesteps.

Using the Jacobian to compute fluxes.

Another advantage of having computed the Jacobian is that, as long as the current state of the atmosphere is sufficiently close to the last state for which we computed \(\mathbb{K}\), Eq. ([jacobian]) can be used to compute the heating rates extremely rapidly. In practice, we use this method as long as \(\mathrm{max}\left(|\hat{T}^{t^{*}}- \hat{T}_\mathbb{K}|\right)\) is below some user defined threshold. Otherwise, the full radiative transfer is computed. This is especially efficient when the atmosphere is approaching equilibrium and requires a large number of small temperature increments to equilibrate the optically thick parts of the atmospheres while the upper atmosphere requires small timesteps to remain stable.

Convergence acceleration

Despite the acceleration procedures described above, the convergence time of an atmosphere can remain prohibitive in some cases. This is usually due to the fact that the upper atmosphere has a short radiative timescale that requires small timesteps for stability while the deep atmosphere is very opaque and evolves very slowly.

To circumvent these issues, the library offers several ways to accelerate the convergence toward an equilibrium state for which \(\hat{H}(\hat{\phi}_\mathrm{eq})=\vec{0}\). It should be stressed, however, that the evolution trajectory followed by the atmosphere toward this state cannot be regarded as a temporal evolution.

For each timestep, we first compute the heating rates for all the processes in the regular way described above over a duration \(\Delta t^{*}\). We then identify radiative zones as groups of adjacent layers where \(H_{\mathrm{n}}^\mathrm{conv} = H_{\mathrm{n}}^\mathrm{madj} = H_{\mathrm{n}}^\mathrm{cond} = H_{\mathrm{n}}^\mathrm{rain} = 0\). A radiative zone can be composed of a single layer. The remaining layers of the atmosphere are grouped in stacks of adjacent layers that we will call convective zones for convenience, even though the energy exchange in these zones might be due to formation and reevaporation of rains.

Acceleration in radiative zones

In radiative zones, we start by defining a base timescale, \(\tau^{*}_\mathrm{b}\), as the smallest radiative timescale of a radiative layer in the atmosphere. Then, for any radiative layer \(\mathrm{n}\), the heating rate is multiplied by \(\tau^{*}_\mathrm{n}/\tau^{*}_\mathrm{b}\) before the timestepping is performed.

In a purely radiative atmosphere this would be equivalent to advancing every layer independently using its own radiative timescale. In practice, this still allows deep radiative zones with very long timescales to converge in as many timesteps as their counterparts in the upper atmosphere.

Acceleration in convective zones

In convective zones, the layers in the zone are not independent as they directly exchange energy in a conservative way: energy that is taken in a part of the zone is redistributed in some other part. Using the scheme above is such zones completely upsets the balance.

To understand why, let us consider a simple 2-layer convective zone where radiation heats the base layer and cools the upper one. When the layer becomes unstable, dry convection carries the surplus of heat from the bottom to the top layer. Equilibrium is reached when the convective energy flux equals both the net radiative heating of the base layers and the net radiative cooling of the top layer. Now, one can see that if the radiative heating rates are multiplied by two independent factors, our equilibrium solution is not an equilibrium anymore and the system will reach a new unphysical equilibrium.

To avoid this problem, we treat an entire convective zone as one layer during the acceleration part. We compute the average radiative heating rate in the zone by computing the net radiative fluxes at the top and bottom of the convective zone and dividing it by its total mass. Then we compute a radiative timescale for the whole zone, \(\tau^{*}_\mathrm{av}\). For the moment, we use the smallest radiative timescale of any layer inside the zone. Finally, each layer receives the same average radiative heating rate multiplied by \(\tau^{*}_\mathrm{av}/\tau^{*}\).

References

[1]

By default, the gravity is kept constant and equal to the surface gravity (\(g_0\)) provided by the user. If a planetary radius is provided, the gravity varies with altitude, \(z\) following \(g=g_0/(1+z/ R_p)^2\).