Modeling atmospheres and their evolution with exo_k: a practical guide

Author: Jeremy Leconte (CNRS/LAB/Univ. Bordeaux)

Since exo_k allows anyone to easily read radiative data and interpolate into them, it seemed silly not to incorporate a module to test this directly in a real-life situation. This is why we implemented the Atm class that will allow you to build atmospheres with any composition, and compute their transmission/emission spectra with any combination of radiative data that the library can handle. It can also compute the radiative heating rates in each atmospheric layer.

But why stop there? Indeed, heating rates allow you to compute the radiative equilibrium of an atmosphere. One thing leading to another, we ended up with a full fledged atmospheric model able to compute the evolution of an atmosphere in time that accounts for radiation (of course), dry convection, but also moist convection, condensation, turbulent diffusion, etc. This model is encapsulated in the :class:~exo_k.atm_evolution.atm_evol.Atm_evolution class that will be presented below.

See also

The goal of this section is to give you a practical tour of how to use these classes in a realistic setting. If you want more information on the inner working of the various modules, you can go to the “Principles of atmospheric modeling with Exo_k” page. For the details on how to load radiative data and some general features about the library, please go see the tutorial on “Getting started with radiative data”.

Initialization

First, let’s make some common initializations. To install the library, see the ‘Getting Exo_k’ page.

Important

For efficiency reasons, the atmospheric models work with a single unit system, which is the international MKS system. To avoid load radiative data in another unit system, we encourage you to always use the global setting xk.Settings().set_mks(True) that will automatically do all the conversions for you.

import exo_k as xk
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
import time,sys,os

xk.Settings().set_mks(True)

We also set up the generic path toward the radiative data

datapath = '../data/'

xk.Settings().set_search_path(datapath + 'xsec', path_type='xtable')
xk.Settings().set_search_path(datapath + 'corrk', path_type='ktable')
xk.Settings().set_search_path(datapath + 'cia', path_type='cia')
xk.Settings().search_path

{'ktable': ['/builds/jleconte/exo_k/data/corrk'],
 'xtable': ['/builds/jleconte/exo_k/data/xsec'],
 'cia': ['/builds/jleconte/exo_k/data/cia'],
 'aerosol': ['/builds/jleconte/exo_k/tutorials']}

# Uncomment the line below if you want to enable interactive plots
#%matplotlib notebook
plt.rcParams["figure.figsize"] = (7,4)
from matplotlib import cycler
colors = cycler('color',[plt.cm.inferno(i) for i in np.linspace(0.1,1,5)])
plt.rc('axes', axisbelow=True, grid=True, labelcolor='dimgray', labelweight='bold', prop_cycle=colors)
plt.rc('grid', linestyle='solid')
plt.rc('xtick', direction='in', color='dimgray')
plt.rc('ytick', direction='in', color='dimgray')
plt.rc('lines', linewidth=1.5)

The Atm class: An atmospheric radiative-transfer model

Whether you want to test various assumptions (will my resolution be enough, are my initial data adequate, what weights should I use, etc.) or directly simulate the transmission/emission spectrum of a planet, we implemented a simple 1d atmosphere model so that you can do just that.

This revolves around the Atm() class, which defines an atmosphere with Pressure, Temperature, and composition profiles. In essence, we are using the structure of the Gas_mix object and add the additional attributes that make an atmosphere more than just some gas with a given pressure, temperature, and composition: a gravity, a planetary radius, etc.

As Xtable() objects can be used in a Kdatabase(), you can do the equivalent of Line by Line radiative transfer.

If you want to start from very high-res spectra, first load them into a Xtable(), and then into a Database.

See also

Check out Atm and its base class Atm_profile for details on arguments and options.

Important

Pressure is the primary coordinate used to build the atmospheric profiles (and not altitude). To keep pressure arrays in increasing order, the atmosphere is therefore described from the top down. The first element of each array describing the atmophere (temperature, composition, etc.) refers to the uppermost layer of the atmosphere. The surface is described by the last element of each array.

Also note that since v1.1.0, there is a new scheme for the computation of atmospheric emission/transmission to improve numerical accuracy. This led to the change of some variable names to instantiate atm objects (Nlev, logplev, and tlev are now Nlay, logplay and tlay).

Loading the radiative data that we will use in our atmospheric model

Molecular absorptions: The Kdatabase

The first thing that we need to do is to create a Kdatabase that will contain the data for all the molecules that we want to include in our atmosphere and for which we have cross sections (whether they are monochromatic cross sections of correlated-k coefficients).

Important

An Atm instance will inherit its spectral range and resolution for the radiative computations from the Kdatabase it is given. So it is important to choose radiative data that cover the wavelength range that you want to model.

As a corrolary, even if you want to model an atmospere for which you have only continuum absorptions (like Rayleigh or Collision Induced Absorptions), you still need to provide a Kdatabase that has at least one molecule in it. If you do not, the model will not be able to know which spectral window to model. We will show later how to turnoff the molecular absorptions.

In the example below we will treat a simple examples with only CO2 and H2O. Other gases, if any, will only contribute through continuum absorptions.

k_db=xk.Kdatabase(['CO2', 'H2O'],'R300_0.3-50mu.ktable.SI')
k_db

The available molecules are:
CO2->/builds/jleconte/exo_k/data/corrk/CO2_R300_0.3-50mu.ktable.SI.h5
H2O->/builds/jleconte/exo_k/data/corrk/H2O_R300_0.3-50mu.ktable.SI.h5
All tables share a common spectral grid
All tables share a common logP-T grid
All tables share a common p unit
All tables share a common kdata unit

Important

The all the tables in the Kdatabase must have the same Wavenumber, Pressure, Temperature (and gauss-point) grids. If this is not the case, refer to the general tutorial to see how to correct that.

Collision induced absorptions

Now we can load CIA data in a CIA_database. If you do not want to include such constributions, do not provide any CIA database when instantiating your Atm object or use cia_database=None.

cia_db=xk.CIAdatabase(molecules=['N2','H2O','H2','He'], mks=True)
cia_db

The available molecule pairs are:
N2-N2->/builds/jleconte/exo_k/data/cia/N2-N2_2011.hdf5
N2-H2->/builds/jleconte/exo_k/data/cia/N2-H2_2011.hdf5
H2O-N2->/builds/jleconte/exo_k/data/cia/H2O-N2_continuumTurbet_MT_CKD3.3.hdf5
H2O-H2O->/builds/jleconte/exo_k/data/cia/H2O-H2O_continuumTurbet_MT_CKD3.3.hdf5
H2-H2->/builds/jleconte/exo_k/data/cia/H2-H2_2011.hdf5
H2-He->/builds/jleconte/exo_k/data/cia/H2-He_2011.hdf5
All tables do NOT have common spectral grid
You will need to run sample before using the database

Just like with the Kdatabase, all the CIA tables must use the same spectral axis (and the same as the Kdatabase) so we usually need to resample our database on the same grid as follows.

cia_db.sample(k_db.wns)
cia_db

The available molecule pairs are:
N2-N2->/builds/jleconte/exo_k/data/cia/N2-N2_2011.hdf5
N2-H2->/builds/jleconte/exo_k/data/cia/N2-H2_2011.hdf5
H2O-N2->/builds/jleconte/exo_k/data/cia/H2O-N2_continuumTurbet_MT_CKD3.3.hdf5
H2O-H2O->/builds/jleconte/exo_k/data/cia/H2O-H2O_continuumTurbet_MT_CKD3.3.hdf5
H2-H2->/builds/jleconte/exo_k/data/cia/H2-H2_2011.hdf5
H2-He->/builds/jleconte/exo_k/data/cia/H2-He_2011.hdf5
All tables share a common spectral grid

Rayleigh scattering

For the moment, the data for Rayleigh scattering are hardcoded inside the library. So there is no need to load them. You’ll just have to use the rayleigh=True keyword when instantiating your Atm object to state that you want to account for rayleigh scattering in any subsequent radiative transfer calculation (even for emission calculations).

This behavior can also be changed afterward with Atm.set_rayleigh(rayleigh=True/False).

Finally, the rayleigh keyword can be used with some methods to locally force the presence/absence of rayleigh scattering.

Building an atmosphere

Simple troposphere/stratosphere parametrization

There are several ways to build an atmosphere. If you just want a simple atmospheric structure with an isothermal stratosphere on top of a dry troposphere with a fixed adiabatic gradient, you can just use the following options

Nlay = 50 # Number of layers
Tsurf = 220. # K => Surface temperature
Tstrat = 100.  # K => Temperature of the stratosphere
psurf = 640. # Pa => Surface pressure
ptop = 1.e0 # Pa => Pressure at the model top
grav = 3.4 # m/s^2 => Surface gravity
rcp = 0.28 # dimless => R/cp
Rp = None # m => Planetary radius
albedo_surf = 0.2 # m => Surface albedo

composition = {'CO2':'background', 'H2O':1.e-6}
# Volume molar concentrations of the various species
# One gas can be set to `background`.
# It's volume mixing ratio will be automatically computed to fill the atmosphere.

atm_mars = xk.Atm(psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             composition = composition, Nlay=Nlay,
             k_database = k_db,
             cia_database = cia_db, # adding a CIA database entails that CIA will be taken into account.
             rayleigh = True, # by default, rayleigh scattering is turned off
                 )

Important

If the sum of the Volume molar concentrations is not equal to one, the rest of the atmosphere is assumed to be transparent and the molar mass of the gas, if needed, is assumed to be equal to the molar mass of the mix you specified (for example, if you specify that your atmosphere is 25% CO2 and 25% H2O, only half of the atmophere will be radiatively active, but the molar mass will still be about 0.031kg/mol).

Warning

The planetary radius is optional. However, if it is not given, the evolution of gravity with height (in transmission and emission) will not be taken into account.

fig,axs=plt.subplots(1,2,sharey=False,figsize=(7,4))
atm_mars.plot_T_profile(axs[0], yscale='log')
atm_mars.plot_T_profile(axs[1], use_altitudes = True)
fig.tight_layout()

If the stratospheric temperature is set to None, the profile will be isothermal at T=Tsurf.

Nlay = 100 # Number of layers
Tsurf = 1000. # K => Surface temperature
Tstrat = None  # K => Temperature of the stratosphere
psurf = 1.e6 # Pa => Surface pressure
ptop = 1.e-4 # Pa => Pressure at the model top
grav = 10. # m/s^2 => Surface gravity
rcp = 0.28 # dimless => R/cp
Rp = 1.*u.Rjup # m => Planetary radius
albedo_surf = 0. # m => Surface albedo

composition = {'H2':'background', 'He':0.09, 'H2O':1.e-3}

atm_hotjup = xk.Atm(psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             composition = composition, Nlay=Nlay,
             k_database = k_db, cia_database = cia_db,
             rayleigh = True,)

fig,axs=plt.subplots(1,2,sharey=False,figsize=(7,4))
atm_hotjup.plot_T_profile(axs[0], yscale='log')
atm_hotjup.plot_T_profile(axs[1], use_altitudes = True)
fig.tight_layout()

Using your own profiles

If you want, you can completely control the profiles you put in your Atm. For this, you just need to provide arrays for logplay, tlay, and the various concentrations.

The atmosphere is defined from the top down so the first value is the top of atmosphere and the last is the surface.

Here is an example for an Earth with a moist troposphere near saturation.

Notice that for the concentrations, you can put either a single value or an array (or a mix), as long as the size of the array is the same as the size of the T-logP arrays.

Nlay = 100 # Number of layers
Tsurf = 300. # K => Surface temperature
Tstrat = 200.  # K => Temperature of the stratosphere
psurf = 1.e5 # Pa => Surface pressure
ptop = 1.e0 # Pa => Pressure at the model top
grav = 9.81 # m/s^2 => Surface gravity
rcp = 0.28 # dimless => R/cp #create the layer presure grid
Rp = 1.*u.Rearth # m => Planetary radius
albedo_surf = 0.2 # m => Surface albedo

logplay=np.linspace(np.log10(ptop),np.log10(psurf), Nlay)
# design your adiabat
tlay=Tsurf*(10**logplay/psurf)**rcp
tlay=np.where(tlay<Tstrat,Tstrat,tlay)

# create your composition. Here, only water is variable and follows saturation.
h2o_thermo_params={'Latent_heat_vaporization': 2.486e6, 'cp_vap': 4.228e3, 'Mvap': 18.e-3,
            'T_ref': 273., 'Psat_ref': 612.}
h2o=xk.Condensing_species(**h2o_thermo_params)
psat_h2o=h2o.Psat(tlay)
x_h2o = psat_h2o / (psat_h2o + 10**logplay)
for i in range(len(x_h2o)-2,-1,-1):
    if x_h2o[i+1] < x_h2o[i]:
        x_h2o[i] = x_h2o[i+1]

composition = {'CO2':4.e-4,'H2O':x_h2o,'N2':'background'}

#Loads the atmosphere and computes some intermediate things
atm_earth=xk.Atm(logplay = logplay, tlay = tlay, grav = grav,
                Rp = Rp, rcp = rcp,
                albedo_surf = albedo_surf,
                composition = composition,
                k_database = k_db, cia_database = cia_db,
                rayleigh = True)

fig,axs=plt.subplots(1,2,sharey=True,figsize=(7,4))
atm_earth.plot_T_profile(axs[0], yscale='log')
axs[1].plot(x_h2o, atm_earth.play)
axs[1].set_xscale('log')
axs[1].set_xlabel('H2O molar concentration')
fig.tight_layout()

Changing attributes of an atmosphere

In an atmosphere, many attributes are inter related and cannot be changed separately. Changing the pressures, for example, will change the mass of the atmosphere, the altitudes, etc.

So, if you want to change some attributes of the instance you are working with after it has been created, it is recommended to use the methods made for that. Here is a non-exhaustive list: set_logPT_profile, set_T_profile(), set_grav(), set_gas(), set_Mgas(), set_rcp(), set_Rp(), set_Rstar(), set_k_database(), set_cia_database(), set_spectral_range(), set_incoming_stellar_flux(), set_internal_flux(), set_surface_albedo(), set_rayleigh().

See the API reference for completeness.

Optional attributes

At the initialisation, a Mgas can be specified to force a molar mass whatever the composition.

With wn_range or wl_range one can limit the wavenumber/wavelength range of future computations.

flux_top_dw (w/m^2) allows you to specify the stellar flux impinging on the top of atmosphere. If you do not specify anything else, the spectral energy distribution will be the one of a blackbody at T=5778K. This temperature can be changed by using the Tstar keyword at the same time. Alternatively, one can provide a Spectrum object with the keyword stellar_spectrum. This spectrum must be in units of power per wavenumber. The actual units are irrelevant as this spectrum will be renormalized with the value of flux_top_dw provided.

Transmission spectroscopy

To compute the transmission spectrum for your atmosphere, you can use the following method. Among the available options, there is rayleigh, Rstar,

Rs=1.*u.Rsun # m => Radius of the star to compute the transit depth

spec=atm_hotjup.transmission_spectrum(Rstar=Rs)

Important

The result of transmission_spectrum() is a Spectrum object onto wich some methods can be used like plotting, binning, and some algeabric operations.

fig,ax = plt.subplots(1,1, sharex=True, sharey=True)
spec.plot_spectrum(ax, xscale='log', label='Native corr-k resolution')
# bin_down_cp makes a deep copy befire binning so that we keep the initial Spectrum instance
spec.bin_down_cp(xk.wavenumber_grid_R(400.,30000.,200)).plot_spectrum(ax, xscale='log', label=r'$R=200$')
spec.bin_down_cp(xk.wavenumber_grid_R(400.,30000.,30)).plot_spectrum(ax, xscale='log', label=r'$R=30$')
ax.set_ylabel('Depth')
ax.legend()
fig.tight_layout()

Notice that the stellar radius is now an attribute of the Atm instance and does not have to be repeated.

Rayleigh scattering can be turned of for any given computation with the rayleigh=False keyword. However, if you want to turn it on/off permanently for the current instance, you can use atm_hotjup.set_rayleigh(rayleigh=True/False).

spec=atm_hotjup.transmission_spectrum(rayleigh = False)

fig,ax = plt.subplots(1,1, sharex=True, sharey=True)
spec.plot_spectrum(ax, xscale='log')
ax.set_ylabel('Depth')
fig.tight_layout()

Emission spectroscopy

Exo_k enables several modes to compute the emission of a planet:

• Atm.emission_spectrum(mu0=, **kwargs): Simple computation using Schwarzschild’s equation with mu0 as the cosine of the observer zenith angle. This mode cannot handle scattering. The output is given as a flux assuming a constant hemispheric intensity, but it can be used to compute the specific intensity at various angles by changin mu0.

• Atm.emission_spectrum(mu_quad_order=, **kwargs): This mode computes the specific intensity at mu_quad_order angles and integrates it to yield the flux. Exact in the non-scattering limit. Cannot handle scattering.

• Atm.emission_spectrum_2stream(**kwargs): This mode computes the upward and downward fluxes using the constant hemispheric mean two-stream approximation of Toon et al. (1989). This mode can handle scattering and provides fluxes at all the layer interfaces (levels), so that radiative heating rates can be computed. In this case, a stellar flux must be specified.

As expected, in a non-scattering medium and with a zero surface albedo, Atm.emission_spectrum and Atm.emission_spectrum_2stream provide the same top of atmosphere outgoing flux almost down to machine precision. The second will be longer as it computes downward fluxes as well.

Computing top of atmosphere fluxes

TOA_flux = atm_mars.emission_spectrum_2stream( \
                wl_range=[0.1,150.0], #can reduce the wavelength range of the computation to save time
                #rayleigh=True, # whether rayleigh scattering is locally accounted for (even for thermal emission))
                    )

Surf_blackbody = atm_mars.surf_bb()
TOA_blackbody = atm_mars.top_bb()
#print("Forward model duration:")
#%timeit -n 2 -r 2 atm_mars.emission_spectrum(integral=True, wl_range=[0.1,150.0])
#%timeit -n 2 -r 2 atm_mars.emission_spectrum_2stream(integral=True, wl_range=[0.1,150.0])

Important

The result of emission_spectrum() is a Spectrum object onto wich some methods can be used like: computing the total flux, plotting, binning, and some algeabric operations (substract of divide spectra).

print("Top of atmosphere flux:", TOA_flux.total, 'W/m^2')
print("Black body at Tsurf over the computed nu range is:",Surf_blackbody.total,'W/m^2')
print('(sigma*Tsurf**4 is:',xk.SIG_SB*atm_mars.tlay[-1]**4,'W/m^2)')

Top of atmosphere flux: 113.29278956850874 W/m^2
Black body at Tsurf over the computed nu range is: 123.82635458875077 W/m^2
(sigma*Tsurf**4 is: 132.8317491952 W/m^2)

fig,axs=plt.subplots(1,2, sharey=False, figsize=(8,4))
Surf_blackbody.plot_spectrum(axs[0], label='Blackbody ($T_{surf}$)')
TOA_blackbody.plot_spectrum(axs[0], label='Blackbody ($T_{strat}$)')
TOA_flux.plot_spectrum(axs[0], label='OLR')
Surf_blackbody.plot_spectrum(axs[1], per_wavenumber=False, label='Blackbody ($T_{surf}$)')
TOA_blackbody.plot_spectrum(axs[1], per_wavenumber=False, label='Blackbody ($T_{strat}$)')
TOA_flux.plot_spectrum(axs[1], per_wavenumber=False, label='OLR')

axs[0].set_ylabel(r'Flux ($W/m^2/cm^{-1}$)')
axs[1].set_ylabel(r'Flux ($W/m^2/\mu m$)')
axs[1].set_xscale('log')
for ax in axs:
    ax.legend()
fig.tight_layout()

By default, the source function is integrated inside each bin. The user can specify integral=False to evaluate the Planck function at the center of the bins only. It is a little faster, but inacurate for large bins. It can be used when using moderate to high-resolution radiative data.

Let’s see the error we would make with this mode (here at a resolution of 300) using algeabric operations on our spectra

TOA_flux_no_integral = atm_mars.emission_spectrum_2stream( \
                integral = False,
                wl_range=[0.1,150.0], #can reduce the wavelength range of the computation to save time
                rayleigh=True, # whether rayleigh scattering is accounted for (even for thermal emission))
                    )
rel_error = ((TOA_flux_no_integral-TOA_flux)/TOA_flux).abs() # .abs() to have the absolute value

fig,ax = plt.subplots(1,1)
rel_error.plot_spectrum(ax, yscale='log')

ax.set_ylabel('Relative error')
fig.tight_layout()

Turning CIAs on and off

# Computes emission flux for the two cases

atm_earth.set_cia_database(None) # remove CIAs
TOA_flux_nocia = atm_earth.emission_spectrum()

atm_earth.set_cia_database(cia_db) # reset CIAs
TOA_flux_cia = atm_earth.emission_spectrum_2stream()

Surf_blackbody = atm_earth.surf_bb()
TOA_blackbody = atm_earth.top_bb()

print("Top of atmosphere flux:")
print("With CIAs   :", TOA_flux_cia.total, 'W/m^2')
print("Without CIAs:", TOA_flux_nocia.total, 'W/m^2')
print("we can see the reduced emission with CIAs")

Top of atmosphere flux:
With CIAs   : 260.76473417977724 W/m^2
Without CIAs: 291.08951048730773 W/m^2
we can see the reduced emission with CIAs

fig,axs=plt.subplots(1,2, sharey=False, figsize=(8,4))
Surf_blackbody.plot_spectrum(axs[0], label='Blackbody ($T_{surf}$)')
TOA_blackbody.plot_spectrum(axs[0], label='Blackbody ($T_{strat}$)')
TOA_flux_nocia.plot_spectrum(axs[0], label=r'OLR w/o CIA')
TOA_flux_cia.plot_spectrum(axs[0], label=r'OLR w/ CIA')

Surf_blackbody.plot_spectrum(axs[1], per_wavenumber=False, label='Blackbody ($T_{surf}$)')
TOA_blackbody.plot_spectrum(axs[1], per_wavenumber=False, label='Blackbody ($T_{strat}$)')
TOA_flux_nocia.plot_spectrum(axs[1], per_wavenumber=False, label=r'OLR w/o CIA')
TOA_flux_cia.plot_spectrum(axs[1], per_wavenumber=False, label=r'OLR w/ CIA')

axs[0].set_ylabel(r'Flux ($W/m^2/cm^{-1}$)')
axs[1].set_ylabel(r'Flux ($W/m^2/\mu m$)')
axs[1].set_xscale('log')
for ax in axs:
    ax.legend()
fig.tight_layout()

Contribution Function

Once an emission spectrum has been computed, one can use the contribution_function() method to compute the spectral contribution function of each layer of the atmosphere. The result is an array of size (Number of radiative layers (Nlay-1), Number of wavenumber bins) that give the contribution function in the bin in W/m^2/str/cm-1.

To integrate over a wavenumber band, one must mutliply by the bin width before adding the values.

In the example below, one can see that in the 15 micron band of CO2, Earth’s atmosphere emision mostly comes from the lower stratsophere whereas the lower troposphere and surface contribute much more in the 10 micron radiative window.

contribution_function = atm_earth.contribution_function()
fig,ax=plt.subplots(1,1,sharey=False,figsize=(3,4))
#for contrib in contribution_function.T:
#    ax.plot(contrib)
ax.plot(contribution_function[:,k_db['CO2'].wlindex(15.)], atm_earth.plev[1:-1], label=r'15 $\mu$m (CO$_2$ band)')
ax.plot(contribution_function[:,k_db['CO2'].wlindex(20.)], atm_earth.plev[1:-1], label=r'10 $\mu$m', linestyle='--')
ax.invert_yaxis()
ax.set_yscale('log')
ax.set_xlabel(r'Contribution Function (W/m$^2$/str/cm$^{-1}$)')
ax.set_ylabel(r'Pressure (Pa)')
ax.legend()
fig.tight_layout()

origin = 'upper'
delta = 0.025

x = k_db.wls
y = atm_earth.play[1:]
X, Y = np.meshgrid(x, y)
Z = atm_earth.exp_minus_tau()

nr, nc = Z.shape

plt.rc('axes',grid=False)
fig, ax = plt.subplots(constrained_layout=True)
ax.invert_yaxis()
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(left=0.3,right=20.)
CS = ax.contourf(X, Y, Z, 10, cmap=plt.cm.bone, origin=origin)
ax.set_title('Atm. Transmittance')
ax.set_xlabel('Wavelength (mu)')
ax.set_ylabel('Log(Pressure/Pa)')
# Make a colorbar for the ContourSet returned by the contourf call.
cbar = fig.colorbar(CS)
cbar.ax.set_ylabel('Exp(-tau)')
plt.rc('axes',grid=True)
# Add the contour line levels to the colorbar

Adding reflected light

To account for reflected light from the star, one just needs to specify the stellar flux impinging on the planet’s top-of-atmosphere (TOA) as well as the spectral energy density of the star.

To do this, there are two different modes

Using a stellar blackbody

For this mode, you just need to specify the bolometric stellar flux at the TOA (in W/m^2) with the flux_top_dw keyword, and the stellar temperature with the Tstar keyword at initilization (or using Atm.set_incoming_stellar_flux()).

Nlay = 50 # Number of layers
Tsurf = 200. # K => Surface temperature
Tstrat = 100.  # K => Temperature of the stratosphere
psurf = 640. # Pa => Surface pressure
ptop = 1.e0 # Pa => Pressure at the model top
grav = 3.4 # m/s^2 => Surface gravity
rcp = 0.28 # dimless => R/cp
Rp = None # m => Planetary radius
albedo_surf = 0.2 # m => Surface albedo
Tstar = 5570. # K => Stellar blackbody temperature
flux_top_dw = 120. # W/m^2 => stellar flux at the TOA.
wn_albedo_cutoff = 2500.

composition = {'CO2':'background', 'H2O':1.e-6}
# Volume molar concentrations of the various species
# One gas can be set to `background`.
# It's volume mixing ratio will be automatically computed to fill the atmosphere.

atm_mars = xk.Atm(psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             wn_albedo_cutoff = wn_albedo_cutoff,
             composition = composition, Nlay=Nlay,
             k_database = k_db, cia_database = cia_db,
             rayleigh = True,
             flux_top_dw = flux_top_dw, Tstar = Tstar)

TOA_flux = atm_mars.emission_spectrum_2stream()

Surf_blackbody = atm_mars.surf_bb()
TOA_blackbody = atm_mars.top_bb()

fig,ax = plt.subplots(1,1, figsize=(5,4))
Surf_blackbody.plot_spectrum(ax, per_wavenumber=False, label='Blackbody ($T_{surf}$)')
TOA_blackbody.plot_spectrum(ax, per_wavenumber=False, label='Blackbody ($T_{strat}$)')
TOA_flux.plot_spectrum(ax, per_wavenumber=False, label='Upward TOA flux')

ax.set_ylabel(r'Flux ($W/m^2/\mu m$)')
ax.set_xscale('log')
ax.legend()
fig.tight_layout()

Important

By default, the surface albedo becomes 0. above 2 microns (5000 cm^-1). You can change that with the wn_albedo_cutoff argument at initialization.

Using an input stellar spectrum

We first have to load a stellar spectrum in a Spectrum object.

For this, we can provide an ascii file with two columns, the first one being the spectral axis, and the second one the radiance in flux per spectral unit.

The code assumes that the spectral quantity (wavelength or wavenumber) and units used are the same for both columns (for example a flux in W/micron if the spectral unit is microns. Fluxes in number of photons are not supported).

The precises units for the flux have no importance as the spectrum will be renormalized when input into the model. So any multiplicative factor will be lost.

stellar_spectrum = xk.Spectrum(filename = '../data/stellar_spectra/solar_spectrum.dat',
                spectral_radiance = True, input_spectral_unit = 'nm')

fig,ax = plt.subplots(1,1, figsize=(5,4))
stellar_spectrum.plot_spectrum(ax)
ax.set_ylabel(r'Flux ($W/m^2/cm^{-1}$)')
ax.set_xscale('log')
fig.tight_layout()

Tstar = None # K => Stellar blackbody temperature

atm_mars = xk.Atm(psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             wn_albedo_cutoff = wn_albedo_cutoff,
             composition = composition, Nlay=Nlay,
             k_database = k_db, cia_database = cia_db,
             rayleigh = True,
             flux_top_dw = flux_top_dw,
             Tstar = Tstar,
             stellar_spectrum = stellar_spectrum,
             )

TOA_flux = atm_mars.emission_spectrum_2stream()

Surf_blackbody = atm_mars.surf_bb()
TOA_blackbody = atm_mars.top_bb()

fig,ax = plt.subplots(1,1, figsize=(5,4))
Surf_blackbody.plot_spectrum(ax, per_wavenumber=False, label='Blackbody ($T_{surf}$)')
TOA_blackbody.plot_spectrum(ax, per_wavenumber=False, label='Blackbody ($T_{strat}$)')
TOA_flux.plot_spectrum(ax, per_wavenumber=False, label='Upward TOA flux')

ax.set_ylabel(r'Flux ($W/m^2/\mu m$)')
ax.set_xscale('log')
ax.legend()
fig.tight_layout()

Heating rates

A by-product of the computation of the various fluxes in the atmosphere is that we know how much radiative energy is deposited and reradiated by each atmospheric layer, i.e. the heating rates.

They can be computed as follows:

heat_rates, net_fluxes = atm_mars.heating_rate()
# Heat_Rate are the heating rates in each layer in W/kg (from the top-down)
# net_fluxes are the net bolometric fluxes (positive upward) in W/m^2 at the top of each layer.

fig,axs=plt.subplots(1,3,sharey=True)
axs[0].set_yscale('log')
axs[0].plot(heat_rates, atm_mars.play, label='Solar Heating')
axs[1].plot(net_fluxes, atm_mars.play, label='Net solar flux')
axs[2].plot(atm_mars.tlay, atm_mars.play)
axs[0].set_xlabel('Heating Rate (W/kg)')
axs[1].set_xlabel(r'Flux (W/m$^2$)')
axs[2].set_xlabel('Temperature (K)')
axs[0].set_ylabel('Pressure (Pa)')
axs[0].invert_yaxis()
fig.subplots_adjust(hspace=0, wspace=0)
fig.tight_layout()

Important

All the energy that is absorbed by the bottom boundary of the model (the surface) is added to the last atmospheric layer (remember that the layers are organized from the top-down), which effectively serves as a surface layer. This is why there often seems to be a discontinuity in the heating rates in the bottom layer.

An internal_flux (Flux coming from the subsurface in W/m^2) can also be specified. In the case, this heating is added to the budget of the surface layer as well.

In our example above, one can easily see from the left panel that the atmosphere is not in radiative equilibrium localy: * the surface receives more energy that it emits, whereas the troposphere below the 1mbar level wants to cool, which should drive convection.

• The heating rates in the stratosphere are positive which means that it is too cold and wants to warm.

Finally, it can be seen in the middle panel that the atmosphere is not is global radiative balance either, a negative net flux meaning that the atmosphere receives more light from its star than it is able to emit.

Bringing the atmosphere back to equilibrium will be the topic of the next section. In the meantime, to have a more complete diagnostic, it is possible to compute the heating rates and fluxes separately for the stellar and self-emission contributions.

Fup_emis, Fdw_emis, Fnet_emis, H_emis,\
 Fup_stel, Fdw_stel, Fnet_stel, H_stel = atm_mars.bolometric_fluxes_2band()

fig,axs=plt.subplots(1,3,sharey=True)
axs[0].set_yscale('log')
axs[0].plot(H_stel, atm_mars.play, label='Stellar')
axs[0].plot(H_emis, atm_mars.play, label='Emis')
axs[1].plot(Fdw_stel, atm_mars.play, label='Stellar down')
axs[1].plot(Fup_stel, atm_mars.play, label='Stellar up')
axs[1].plot(Fdw_emis, atm_mars.play, label='Emis down')
axs[1].plot(Fup_emis, atm_mars.play, label='Emis up')
axs[2].plot(Fnet_stel, atm_mars.play, label='Stellar net')
axs[2].plot(Fnet_emis, atm_mars.play, label='Emis net')
axs[0].set_xlabel('Heating Rate (W/kg)')
axs[1].set_xlabel(r'Flux (W/m$^2$)')
axs[2].set_xlabel('Temperature (K)')
axs[0].set_ylabel('Pressure (Pa)')
axs[0].legend(loc='upper left', fontsize='small')
axs[1].legend(loc='best', fontsize='small')
axs[2].legend(loc='best', fontsize='small')
axs[0].set_xlim((-0.3,0.3))
axs[0].invert_yaxis()
fig.tight_layout()

Important

In our conventions, upward (up) and net fluxes are counted positive upward and downward (dw) fluxes are positive downward. Hence, Fnet = Fup - Fdw.

The total heating rates and fluxes are the sum of the contribution of the stellar and self-emission contributions.

Modeling atmospheric evolution with Atm_evolution

Now that we know how to model the radiative transfer (and heating rates) for a given atmopspheric structure and composition, it is easy to predict how it evolves in time.

In Exo_k, this is handled by the Atm_evolution class that we present in this section.

Tailoring radiative data

During an evolution, the radiative transfer may be called hundreds of thousands of times. Furthermore, we do not need a high resolution to compute precise heating rates.

Therefore, to gain computation time, we will reduce the resolution of our Ktables. We will also make sure to extend its spectral range to capture the whole planetary and stellar spectral range.

For the demonstration, we will use a really coarse resolution, but one should of course test how the spectral grid affect their results to choose an appropriate setup.

Res = 8.#.4 # Using very low resolution for sake of speed
wn_grid_LR = xk.wavenumber_grid_R(100.,35000.,Res)
print(wn_grid_LR)

k_db = xk.Kdatabase(['CO2', 'H2O'],'R300_0.3-50mu.ktable.SI')
k_db_LR = k_db.bin_down_cp(wnedges = wn_grid_LR, remove_zeros = True)

cia_db_LR=xk.CIAdatabase(molecules=['N2','H2O','H2','He'], mks=True)
cia_db_LR.sample(k_db_LR.wns)

[  100.           113.31484531   128.40254167   145.49914146
   164.87212707   186.82459574   211.70000166   239.8875294
   271.82818285   308.02168489   349.03429575   395.50767229
   448.16890703   507.84190372   575.4602676    652.08191203
   738.90560989   837.28974881   948.77358364  1075.10131861
  1218.24939607  1380.45741861  1564.26318842  1772.54241215
  2008.55369232  2275.98950935  2579.03399172  2922.42837812
  3311.54519587  3752.47231596  4252.10820001  4818.26982911
  5459.81500331  6186.78092504  7010.54123467  7943.98395523
  9001.71313005 10200.27730827 11558.42845272 13097.41532108
 14841.31591026 16817.41416518 19056.62684586 21593.98723106
 24469.19322642 27727.22845231 31419.06602857 35000.        ]

fig,ax = plt.subplots(1,1, figsize=(5,4))
k_db['H2O'].plot_spectrum(ax, p=1., t=300., g=1., xscale = 'log', yscale = 'log')
k_db_LR['H2O'].plot_spectrum(ax, p=1., t=300., g=1., xscale = 'log', yscale = 'log')
ax.set_ylim((1.e-34, 1.e-20))
ax.set_xlim((0.4, 10.))
fig.tight_layout()

Initializing an atmospheric evolution

Except for a few exceptions that we will discuss in more details, the initialization of an Atm_evolution instance requires the same arguments as the initialization of an Atm object. The initial temperature profile can be specified exactly the same way, with the same options.

The most notable exceptions are: * bg_vmr: The molar concentrations (sometimes called the volume mixing ratios) of the various molecules constituting the background gas. This replaces the composition argument in a standard Atm object, but takes the same type of input. The reason for this change of name will be more evident in the next sections when we include active tracers (such as condensible species). bg_vmr can receive a dictionary of arrays if one wants to model an atmosphere with a fixed (but not uniform) distribution of gases. Variable compositions will be treated later. * cp: Specific heat capacity of the gas at constant pressure (J/kg/K), that is needed to compute a temporal evolution. If it is not given, the code will try to compute it for some known molecules at a standard temperature. * M_bg: molar mass of the background gas (kg/mol; optional). It is possible to use a fixed molar mass of the background gas instead of the one computed from the bg_vmr.

See also

Check out Atm_evolution for details on arguments and options.

Nlay = 50 # Number of layers
Tsurf = 200. # K => Surface temperature
Tstrat = 100.  # K => Temperature of the stratosphere
psurf = 640. # Pa => Surface pressure
ptop = 1.e0 # Pa => Pressure at the model top
grav = 3.4 # m/s^2 => Surface gravity
rcp = 0.28 # dimless => R/cp
Rp = None # m => Planetary radius
albedo_surf = 0.2 # m => Surface albedo
Tstar = 5570. # K => Stellar blackbody temperature
flux_top_dw = 120. # W/m^2 => stellar flux at the TOA.

composition = {'CO2':'background', 'H2O':1.e-6}
# Volume molar concentrations of the various species
# One gas can be set to `background`.
# It's volume mixing ratio will be automatically computed to fill the atmosphere.

evol_mars = xk.Atm_evolution(Nlay=Nlay, psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             bg_vmr = composition, # notice the change in argument name
             k_database = k_db_LR, cia_database = cia_db_LR,
             rayleigh = True,
             flux_top_dw = flux_top_dw, Tstar = Tstar)

Important

An Atm_evolution instance always as the current atmospheric structure as an Atm instance that is accessible through the .atm attribute so that we can use all the methods that we have seen above to diagnose the current state of the atmosphere, such as saving and plotting the initial TP profile.

t_ini = evol_mars.atm.tlay
p_ini = evol_mars.atm.play

fig,ax=plt.subplots(1,1,figsize=(3,4))
evol_mars.atm.plot_T_profile(ax, yscale='log')
fig.tight_layout()

Setting evolution options

Now, that we have an initial atmospheric structure, we must decide what processes we will take into account when modeling its evolution, along with some parameters.

Here is a (non-exhaustive) list of the (frequently evolving) available options: * rayleigh (bool, default: True) * convection (bool, default: False) * convective_transport (bool, default: True) * diffusion (bool, default: False) * condensation (bool, default: False) * rain (bool, default: False) * latent_heating (bool, default: True) * moist_convection (bool, default: False) * moist_inhibition (bool, default: False) * dTmax_use_kernel (float, default: 10.): Maximum temperature difference above which the full radiative transfer is recomputed. A small value increases the computation time. Too large a value may result in numerical instabilities. (See the theoretical description of the evolution model). * acceleration_mode (int, default: 0): What kind of physical acceleration scheme you want to use. 0 means no acceleration. 4 is the most aggressive.

These options can be set at initialization or changed at any time through the ``set_options`` method.

evol_mars.set_options(convection=True, dTmax_use_kernel = 0.,
                flux_top_dw = 130.,
                internal_flux = 10.,)

The set_options method can also be used to change some attributes of the underlying atmosphere and/or radiative transfer, like the radiative data used, or whether one wants to include the CIA, the stellar spectrum, etc. As these changes entail the re-initialization of the the underlying Atm instance, including the recomputation of the radiative kernel, it is more efficient to do all these changes in a single call to set_options or even to set them directly at the initialization of the Atm_evolution object.

Computing a time evolution

To move the model forward in time, you just need to call the evolve(N_timestep=) method where N_timestep is the number of timesteps that you want to run.

The duration of those timesteps will be decided by the model to be the largest stable timestep. If the user wants to reduce the timestep, it can specify a timestep_factor which must be <=1. A larger timestep_factor usually (although not always) results in unstable evolutions.

You can call evolve() as many times as you want. It will always pick-up just where you left.

Warning

Due to various optimizations, running twice evolve() for N timesteps or once for 2N timesteps will not (in general) yield exactly the same results. If you want your results to be reproducible to machine precision, be sure to proceed exactly the same way.

def plot_evolution(model, figsize = (6,6), write = False):
    """A small routine to make standard plots of the evolution.
    It returns the figure and axis obsjects.
    Customize at will.
    """
    time_evo = model.time_hist/(xk.DAY)
    if write :
        print('evol_time (days):', time_evo[-1])
        print('last timestep:', (time_evo[-1]-time_evo[-2]),'day')
        print('Number of steps since last kernel computation:', model.N_last_ker)
    fig,axs=plt.subplots(2,1, sharex=True, figsize = figsize)
    axs[0].plot(time_evo, model.Fnet_top, label='Top of Atmosphere')
    axs[0].plot(time_evo, np.ones_like(model.Fnet_top)*model.atm.internal_flux,
                label='internal flux')
    axs[1].plot(time_evo, model.tlay_hist[:,-1], label='internal flux')
    axs[0].set_ylabel(r'Net Flux')
    axs[1].set_ylabel(r'$T_{surf}$')
    axs[1].set_xlabel(r'Time (days)')
    axs[0].legend()
    fig.tight_layout()
    return fig, axs

model = evol_mars
start=time.time()

model.evolve(N_timestep=100) # where the magic happens

print('simulation time:',time.time()-start,'s')

fig, axs = plot_evolution(model, write = True)

simulation time: 4.156727075576782 s
evol_time (days): 5.910939175052595
last timestep: 0.04653111662635823 day
Number of steps since last kernel computation: 100

One can see that the surface is warming up as the Net TOA flux is smaller than the internal flux (in this case, smaller than 0.).

Plotting profiles

We can also show the profiles for the: * Heating rates: in attribute H_ave[process] where process is a string that can be one of ‘rad’ (radiation), ‘conv’ (dry convection), ‘cond’ (condensation), ‘madj’ (moist convection), ‘rain’ (rain evaporation), ‘tot’ (total). They are given in W/kg and can be converted in K/s by dividing by self.cp. * Net fluxes: in attribute Fnet[process] where process is as above. They are given at the top of each layer in W/m^2. * Temperature and pressure through self.atm.tlay and self.atm.play

def plot_profiles(model, tracers = False, figsize=(7,4),
                  processes=None, mol = 'H2O'):
    """A small routine to make standard plots of the atmospheric profiles.
    It returns the figure and axis obsjects.
    Customize at will.
    """
    if tracers:
        fig,axs=plt.subplots(1,4, sharey=True, figsize=figsize)
        if processes is None:
            processes = ['rad','conv','madj','cond','rain','tot']
    else:
        fig,axs=plt.subplots(1,3, sharey=True, figsize=figsize)
        if processes is None:
            processes = ['rad','conv','tot']
    for label in processes:
        axs[0].plot(model.heating_rate(label)*xk.DAY/model.cp, model.atm.play, label=label)
        axs[1].plot(model.net_flux(label), model.atm.play, label=label)
    axs[2].plot(model.atm.tlay, model.atm.play, label='Final')
    axs[0].set_ylabel(r'Pressure')
    axs[0].set_xlabel(r'Heating (K/day)')
    axs[1].set_xlabel(r'Fluxes')
    axs[2].set_xlabel(r'T(K)')
    if tracers:
        axs[3].plot(model.tracers[mol], model.atm.play, label=mol)
        axs[3].plot(model.qsat(mol), model.atm.play, label=r'$q_{sat}$', linestyle='--')
        axs[3].set_xlabel(r'q (kg/kg)')
        axs[3].set_xscale('log')
    for ax in axs:
        ax.legend(fontsize='small', loc='upper right')
        ax.set_yscale('log')
    axs[0].invert_yaxis()
    fig.tight_layout()
    return fig, axs

def plot_profiles(model, fig=None, axs=None, tracers = False, rh=False, figsize=(7,4),
                  processes=None, mol='H2O', label_suffix=''):
    """A small routine to make standard plots of the atmospheric profiles.
    It returns the figure and axis obsjects.
    Customize at will.
    """
    if (fig is None) and (axs is None):
        if tracers:
            fig_loc, axs_loc = plt.subplots(1,4, sharey=True, figsize=figsize)
            if processes is None:
                processes = ['rad','conv','madj','cond','rain','tot']
        else:
            fig_loc,axs_loc=plt.subplots(1,3, sharey=True, figsize=figsize)
            if processes is None:
                processes = ['rad','conv','tot']
    else:
        fig_loc, axs_loc = fig, axs
    for label in processes:
        axs_loc[0].plot(model.heating_rate(label)*xk.DAY/model.cp, model.atm.play, label=label+label_suffix)
        axs_loc[1].plot(model.net_flux(label), model.atm.play, label=label+label_suffix)
    axs_loc[2].plot(model.atm.tlay, model.atm.play, label='Final'+label_suffix)
    axs_loc[0].set_ylabel(r'Pressure')
    axs_loc[0].set_xlabel(r'Heating (K/day)')
    axs_loc[1].set_xlabel(r'Fluxes')
    axs_loc[2].set_xlabel(r'T(K)')
    if tracers:
        if rh:
            axs_loc[3].plot(model.tracers[mol]/model.qsat(mol), model.atm.play, label=mol+label_suffix)
            axs_loc[3].set_xlabel(r'RH')
        else:
            axs_loc[3].plot(model.tracers[mol], model.atm.play, label=mol+label_suffix)
            axs_loc[3].plot(model.qsat(mol), model.atm.play, label=r'$q_{sat}$'+label_suffix, linestyle='--')
            axs_loc[3].set_xlabel(r'q (kg/kg)')
            axs_loc[3].set_xscale('log')
    for ax in axs_loc:
        ax.legend(fontsize='small', loc='upper right')
        ax.set_yscale('log')
    if fig is None:
        axs_loc[0].invert_yaxis()
        fig_loc.tight_layout()
    return fig_loc, axs_loc

fig, axs = plot_profiles(evol_mars)
axs[2].plot(t_ini, p_ini, label='Ini')
_ = axs[2].legend(fontsize='small', loc='upper right')

As anticipated, the stratosphere was too cold and is now getting warmer.

If evolve is called again, the evolution will resume from its last point so that you can iteratively reach equilibrium.

Acceleration

Using the radiative kernel

Most of the time, we are only concerned with the equilibrium state of the atmosphere. However, it can take lots of iterations to get there. The evolution model thus includes several ways of accelerating the computation.

The easiest of all is to allow the model to use a linear approximation of a past radiative transfer computation if the temperature profile has not changed by more than dTmax_use_kernel Kelvin.

This can be changed as follows.

evol_mars.set_options(dTmax_use_kernel = 10.)

model = evol_mars
start=time.time()
model.evolve(N_timestep=1000)

print('simulation time:',time.time()-start,'s')
fig, axs = plot_evolution(model, write = True)

simulation time: 2.3474905490875244 s
evol_time (days): 46.50681636282888
last timestep: 0.046468272314868386 day
Number of steps since last kernel computation: 148

Notice that the simulation ran for almost as much time as before, but we have simulated a period about ten times as long (10 times the number of iterations).

The price to pay are those sawtooth patterns that appear each time that we compute the kernel because our linear approximation was not that good anymore. These however do not affect the final equilibrium state as long as the model converges.

If dTmax_use_kernel is too large, the model may not converge.

Other acceleration modes

To reach equilibrium even faster it is possible to use other acceleration modes. They are described in the ‘Principles of atmospheric modeling with Exo_k’ page. They can be activated through the acceleration_mode argument as follows.

The most efficient one to date is acceleration_mode = 4 which accelerates both convective and radiative zones. Use acceleration_mode = 3 if you only want to accelerate radiative zones.

In this case and when we are close to equilibrium, it might be even faster to force the recomputation of the kernel a few times during an evolution step. This can be done through the evolve(N_kernel = ) argument that specifies a maximum number of step between two computations of the kernel.

To show how fast this is, let’s compute the equilibrium for a new stellar flux.

evol_mars.set_options(flux_top_dw = 200., dTmax_use_kernel = 10., acceleration_mode = 4)

model = evol_mars
start=time.time()

model.evolve(N_timestep = 40, N_kernel = 10)

print('simulation time:',time.time()-start,'s')
fig, axs = plot_evolution(model, write = True)

simulation time: 3.4685513973236084 s
evol_time (days): 1.4575971498506874
last timestep: 0.034215043606686146 day
Number of steps since last kernel computation: 36

In only 40 iterations, we are almost in thermal equilibrium!

Saving and loading an atmospheric model

An atmospheric model can be saved to disc using serialization with:

evol_mars.write_pickle('evol_mars.pickle')

Now, one can immediately load it again. It will have all the of its attributes and options saved so that one can directly see the current state as well as make it evolve again.

evol_mars2 = xk.load_atm_evolution('evol_mars.pickle')

model = evol_mars2
try:
    model.evolve(N_timestep = 40, N_kernel = 10)
    print('I have all I need to evolve again')
except:
    print('This Atm_evolution instance is not linked to any k_database')
    print('We will reset the radiative model before running again')
    model.set_options(k_database = k_db_LR, cia_database = cia_db_LR)
    model.evolve(N_timestep = 40, N_kernel = 10)
print('simulation time:',time.time()-start,'s')
fig, axs = plot_evolution(model, write = True)

I have all I need to evolve again
simulation time: 6.15186071395874 s
evol_time (days): 1.3682041583566824
last timestep: 0.03419611852372939 day
Number of steps since last kernel computation: 40

fig, axs = plot_profiles(evol_mars2)
axs[2].plot(t_ini, p_ini, label='Ini')
_ = axs[2].legend(fontsize='small', loc='upper right')

However, when saving it this way, the model may take some space on disk as it saves along all the opacity databases.

If one wants to save only the current state and the options, it is possible to use the data_reduction_level=2 argument in write_pickle.

However, if one wants to make the atmosphere evolve again, they will need to reset the radiative model with the opacities as follows

evol_mars.write_pickle('evol_mars_small.pickle', data_reduction_level=2)

evol_mars3 = xk.load_atm_evolution('evol_mars_small.pickle')

model = evol_mars3
try:
    model.evolve(N_timestep = 40, N_kernel = 10)
    print('I have all I need to evolve again')
except:
    print('We will reset the radiative model before running again')
    model.set_options(k_database = k_db_LR, cia_database = cia_db_LR)
    model.evolve(N_timestep = 40, N_kernel = 10)
print('simulation time:',time.time()-start,'s')
fig, axs = plot_evolution(model, write = True)

This Atm_evolution instance is not linked to any k_database
Use self.set_options(k_database=, ...)
We will reset the radiative model before running again
simulation time: 9.829702854156494 s
evol_time (days): 1.368186821276735
last timestep: 0.03419611689270208 day
Number of steps since last kernel computation: 40

Equilibrating an atmosphere

While you can run evolve as many times as you want to reach equilibrium, the library also includes an equilibrate method that directly runs evolve until thermal equilibrium is met.

This method uses the following arguments: * Fnet_tolerance (W/m^2): The model will be considered equilibrated when the net TOA flux (corrected for the internal flux) will remain between + and - Fnet_tolerance during a call to evolve. * N_iter_max: the maximum number of calls to evolve. * N_timestep_ini: The number of timesteps in the first call to evolve. The number of timesteps will then be doubled each time until it reaches N_timestep_max. * verbose: to show steps and diagnostics. * Any keyword argument for evolve can be passed.

Some atmospheres have an intrinsic variability. Fnet_tolerance must be large enough to account for this variability.

After a call to equilibrate, only the diagnostics for the last call to evolve are available (e.g. heating rates, fluxes, etc.)

Important

For the moment, there is no convergence criterion involving the tracer distribution. Hence, after equilibrate, a model may be in thermal equilibrium but the tracer distribution might not be stationary.

Nlay = 50 # Number of layers
Tsurf = 200. # K => Surface temperature
Tstrat = 100.  # K => Temperature of the stratosphere
psurf = 640. # Pa => Surface pressure
ptop = 1.e0 # Pa => Pressure at the model top
grav = 3.4 # m/s^2 => Surface gravity
rcp = 0.28 # dimless => R/cp
Rp = None # m => Planetary radius
albedo_surf = 0.2 # m => Surface albedo
Tstar = 5570. # K => Stellar blackbody temperature
flux_top_dw = 120. # W/m^2 => stellar flux at the TOA.

composition = {'CO2':'background', 'H2O':1.e-6}
# Volume molar concentrations of the various species
# One gas can be set to `background`.
# It's volume mixing ratio will be automatically computed to fill the atmosphere.

evol_mars = xk.Atm_evolution(Nlay=Nlay, psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             bg_vmr = composition, # notice the change in argument name
             k_database = k_db_LR, cia_database = cia_db_LR,
             rayleigh = True,
             flux_top_dw = flux_top_dw, Tstar = Tstar,
             convection = True,
            )

evol_mars.equilibrate(Fnet_tolerance = 0.1, verbose = True)
fig, axs = plot_evolution(evol_mars, write = True)
fig, axs = plot_profiles(evol_mars)

iter: 1, N_timestep: 100
Fnet mean: -2.07 W/m^2, (min:-13.8, max:-0.544)
timestep: 0.0531 d | 4.59e+03 s, total time: 0.0148 yr
iter: 2, N_timestep: 200
Fnet mean: -0.476 W/m^2, (min:-0.552, max:-0.405)
timestep: 0.0531 d | 4.58e+03 s, total time: 0.0291 yr
iter: 3, N_timestep: 400
Fnet mean: -0.28 W/m^2, (min:-0.434, max:-0.161)
timestep: 0.053 d | 4.58e+03 s, total time: 0.0581 yr
iter: 4, N_timestep: 800
Fnet mean: -0.0612 W/m^2, (min:-0.201, max:-0.0111)
timestep: 0.053 d | 4.58e+03 s, total time: 0.116 yr
iter: 5, N_timestep: 1600
Fnet mean: -0.00405 W/m^2, (min:-0.0632, max:-2.68e-05)
timestep: 0.053 d | 4.58e+03 s, total time: 0.232 yr
evol_time (days): 84.78933052052984
last timestep: 0.052972794704771786 day
Number of steps since last kernel computation: 1600

Tracers and condensing species

Putting tracers in your model

Now let’s try to model an Earth-like planet with some water vapor in the atmosphere.

The first step is to add one or more tracers (at least two if you want a condensing species) through the tracers argument. The input is a dictionary whose keys are the tracer names and the values are another dictionary containing the properties of the tracer. In particular, for a condensing species, you want to specify that it is a 'type':'vapor' and the name of the condensed tracer with 'condensed_form':'H2O_liq' (This tracer must also be added. However, it does not need any specific property).

There are several ways to specify the initial values/surface sources for the various tracers, as will be shown below. For the moment, let’s just add a constant initial mass concentration (kg/kg of gas) profile through 'init_value': q_h2o.

As we want the code to be as generic as possible, we must now specify the thermodynamical properties of the condensing pairs we want to model through the condensing_species argument. There is an example for water vapor below.

Thanks to this mechanism, there is no limitation on the number of condensing species that one can add, as long as we know the thermodynamical properties.

Important

The name of the vapor tracer must be the same as the name of the condensing species. This is how the library links the two. Moreover, if you want the vapor to be radiatively active, the name must be the same as the one in the Kdatabase. For example, here, putting ‘H2O_vap’ in tracers and condensible species would work, but water vapor would have no effect on the radiative transfer.

Warning

The tracers are parametrized in terms of their specific concentration, q, meaning their mass per mass of gas (in kg/kg), not their molar concentration or volume mixing ratios. This should not be mistaken with the mass mixing ratio that is the mass of vapor per mass of dry air. It should be noted that the specific concentration is often called mass mixing ratio in many contexts.

# a simple function to test water conservation
def water_content(model):
    idx_vap, idx_cond = model.condensing_pairs_idx[model.condensing_species_idx['H2O']]
    return np.sum(model.tracers.qarray[idx_vap]*model.atm.dmass) \
      +np.sum(model.tracers.qarray[idx_cond]*model.atm.dmass) \
      +model.tracers.qsurf[idx_cond]

Nlay = 50 # Number of layers
Tsurf = 280. # K => Surface temperature
Tstrat = 150.  # K => Temperature of the stratosphere
psurf = 1.e5 # Pa => Surface pressure
ptop = 10.e0 # Pa => Pressure at the model top
grav = 9.81 # m/s^2 => Surface gravity
rcp = 0.28 # dimless => R/cp
Rp = None # m => Planetary radius
albedo_surf = 0.35 # m => Surface albedo
Tstar = 5570. # K => Stellar blackbody temperature
flux_top_dw = 340. # W/m^2 => stellar flux at the TOA.

composition = {'N2': 'background', 'CO2':340.e-6}

qH2O = 0.0001 # kg/kg of gas

# the tracers dictionary contains the tracers with their names and attributes.
# Here we have two: one for water vapor and one for the condensed phase.
# !!! The vapor one needs to have the same name as in the k_database if
#      the user wants to include its radiative effect !!!
tracers = {'H2O':{'type':'vapor', # a 'vapor' type specifies that it can condense
                  'condensed_form':'H2O_liq', # if 'type'='vapor' then we need to specify the condensed form
                  'init_value':qH2O, # optional, the initial specific concentration (in kg/kg of total gas)
                 },
         'H2O_liq':{}, # The tracer for the condensed form does not need any attribute
         }

# The condensing_species dictionary captures the necessary thermodynamical
# constant to compute the saturation curves and latent heat.
# There needs to be as many condensing_species as there are tracers with a vapor type.
condensing_species={'H2O':{'Latent_heat_vaporization': 2.086e6, # Latent heat at ref temperature (J/kg)
            'cp_vap': 4.228e3, # Specific heat of the vapor (J/kg/K)
            'Mvap': 18.e-3, # Molar mass of the vapor
            'T_ref': 273., # Reference temperature
            'Psat_ref': 612., # Saturation pressure at reference temeprature
                          }}

tracer_values = None # Dict with arrays of values for the various tracers.
                 # It is not needed here as we start with a constant value
                 # as specified above

evol_earth = xk.Atm_evolution(Nlay=Nlay, psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             bg_vmr = composition, # notice the change in argument name
             k_database = k_db_LR, cia_database = cia_db_LR,
             rayleigh = True,
             flux_top_dw = flux_top_dw, Tstar = Tstar,
             tracers = tracers, condensing_species = condensing_species,
             tracer_values = tracer_values)


initial_water = water_content(evol_earth)
print("water content = ", initial_water, "kg/m^2")

water content =  1.0192660550458716 kg/m^2

Tracer related options

When (condensing) tracers are added, there are several additional options you can set. In particular: * Know that convection=True will insantaneously mix all tracers in convective zones. * condensation=True will condense vapors to bring any layer as close as possible to saturation equilibrium. * moist_convection=True calls a moist-convection adjustement scheme. * The effect of mean molecular weight on both dry and moist convection can be turned on with moist_inhibition=True. * rain=True collects all the condensates left in the atmosphere by the previous processes and carries them downward until they are either re-evaporated or they reach the surface. * diffusion=True applies turbulent diffusive mixing to the tracers using a constant eddy mixing coefficient that is specified using Kzz=value (in m^2/s). This diffusion has no effect on temperature.

# The moist processes are turned off by default and need to be turned on
evol_earth.set_options(convection = True, condensation = True,
            moist_convection = True, moist_inhibition = True,
            diffusion = True, Kzz=1.e-2, # in m^2/s
            rain = True,
            acceleration_mode=0,
            )

Important

When moist processes are included, the physical timescale for those processes can be very different from the radiative timescale, which is the one used to compute the timestep. If the moist processes’ timescale is shorter, the user needs to reduce the timestep with the timestep_factor argument in Atm_evolution.evolve(). The actual timestep computed by the code will be multiplied by timestep_factor before the step is taken. Using timestep_factor>1 may result in unstable evolutions.

model = evol_earth
start=time.time()
model.evolve(N_timestep=4000, timestep_factor=.01)
print('simulation time:',time.time()-start,'s')
print("water content:", water_content(model), "kg/m^2")
print("relative water mass loss:", (water_content(model)-initial_water)/initial_water)
fig, axs = plot_evolution(model, write = True)

simulation time: 17.37640118598938 s
water content: 1.019266055046121 kg/m^2
relative water mass loss: 2.446427898745122e-13
evol_time (days): 112.79747220225067
last timestep: 0.03010573793422111 day
Number of steps since last kernel computation: 3045

fig, axs = plot_profiles(evol_earth, tracers = True)

Specifying initial tracer profiles and tracer sources

For the moment, there are two ways to specify the initial tracer distribution: * As seen above, 'init_value' in the tracers dictionary initializes the tracer with a constant specific concentration profile. * A vertical profile of specific concentration can be provided through tracer_values={'name_of_tracer: array}. The array must have right size and be ordered from the top down, just like the other arrays.

If nothing else is specified, the model will evolve assuming that there is no surface source for tracers.

Tracer sources at the bottom boundary can be modelled with two different approaches: * By putting 'q_deep': value as a property of a given vapor tracer in tracers, you force the specific concentration for this vapor to remain fixed at this value in the surface layer. This is particularly usefull to model gaseous planets with no definite surface. * By putting 'surface_reservoir': value as a property of a given condesate tracer in tracers, you specify that there is a reservoir containg value kg/m^2 of the condensed species at the surface. The atmosphere can freely condense/evaporate mater to and from this reservoir depending on the thermodynamical conditions. This is particularly usefull to model a surface ocean.

Setting up an ocean

The only line that changes to put an ocean is 'H2O_liq':{'surface_reservoir': 1000.},. Here we limit it to 1m (1000. kg/m^2).

tracers = {'H2O':{'type':'vapor', # a 'vapor' type specifies that it can condense
                  'condensed_form':'H2O_liq', # if 'type'='vapor' then we need to specify the condensed form
                  'init_value':0., # we start without any atm vapor
                 },
         'H2O_liq':{'surface_reservoir': 1000.}, # here we put a 1m deep ocean
         }

evol_earth = xk.Atm_evolution(Nlay=Nlay, psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             bg_vmr = composition, # notice the change in argument name
             k_database = k_db_LR, cia_database = cia_db_LR,
             rayleigh = True,
             flux_top_dw = flux_top_dw, Tstar = Tstar,
             tracers = tracers, condensing_species = condensing_species,
             tracer_values = tracer_values,
             # evolution options
             convection = True, condensation = True,
             moist_convection = True, moist_inhibition = True,
             diffusion = True, Kzz = 1.e-1, # in m^2/s
             rain = True,
             acceleration_mode = 0,
             mass_redistribution = False)
initial_water = water_content(evol_earth)
print("water content = ", initial_water, "kg/m^2")

water content =  1000.0 kg/m^2

model = evol_earth
start=time.time()
model.evolve(N_timestep=2000, N_kernel = 1000, timestep_factor=0.9)
print('simulation time:',time.time()-start,'s')
print("water content:", water_content(model), "kg/m^2")
print("relative water mass loss:", (water_content(model)-initial_water)/initial_water)
fig, axs = plot_evolution(model, write = True)
fig, axs = plot_profiles(evol_earth, tracers = True)

simulation time: 2.3483967781066895 s
water content: 999.9999999999677 kg/m^2
relative water mass loss: -3.228706191293895e-14
evol_time (days): 7114.40207750129
last timestep: 3.567518697239393 day
Number of steps since last kernel computation: 1957

Forcing surface concentration

Here, the important option is 'q_deep':1.e-3.

Of course, in this setup, the tracer mass is not conserved as we have an open system.

Also, in this setting, you better make sure that the tracer does not condense at the surface or the model will be unstable.

tracers = {'H2O':{'type':'vapor', # a 'vapor' type specifies that it can condense
                  'condensed_form':'H2O_liq', # if 'type'='vapor' then we need to specify the condensed form
                  'init_value':0., # optional, the initial specific concentration (in kg/kg of total gas)
                  'q_deep':1.e-3,
                 },
         'H2O_liq':{}, # The tracer for the condensed form does not need any attribute
         }

evol_earth = xk.Atm_evolution(Nlay=Nlay, psurf=psurf, ptop=ptop,
             Tsurf=Tsurf, Tstrat=Tstrat,
             grav=grav, rcp = rcp, Rp=Rp,
             albedo_surf = albedo_surf,
             bg_vmr = composition, # notice the change in argument name
             k_database = k_db_LR, cia_database = cia_db_LR,
             rayleigh = True,
             flux_top_dw = flux_top_dw, Tstar = Tstar,
             tracers = tracers, condensing_species = condensing_species,
             tracer_values = tracer_values,
             # evolution options
             convection = True, condensation = True,
             moist_convection = True, moist_inhibition = True,
             diffusion = True, Kzz = 1.e1, # in m^2/s
             rain = True,
             acceleration_mode = 0,
             mass_redistribution = False)
initial_water = water_content(evol_earth)
print("water content = ", initial_water, "kg/m^2")

water content =  0.0 kg/m^2

model = evol_earth
start=time.time()
model.evolve(N_timestep=2000, N_kernel = 1000, timestep_factor=0.01)
print('simulation time:',time.time()-start,'s')
print("water content:", water_content(model), "kg/m^2")
fig, axs = plot_evolution(model, write = True)
fig, axs = plot_profiles(evol_earth, tracers = True, rh = False)

simulation time: 1.7908625602722168 s
water content: 5.002828324240059 kg/m^2
evol_time (days): 67.88828125835565
last timestep: 0.03324392465508197 day
Number of steps since last kernel computation: 1889