Deep-Jupiter TP adiabat with KCl and MgSiO3
This runnable example constructs a one-dimensional, hydrostatic, reversible adiabat for an illustrative deep-Jupiter mixture. It includes explicit H2 and He, KCl condensation, and Mg + SiH4 + 3H2O <=> MgSiO3(s) + 5H2. Both condensates remain suspended; there is no rainout, mixing, radiative transfer or external heat source.
The YAML defines the chemistry and column boundary conditions. The Python
script calls ThermoX.forward for every temperature trial, then finds the
temperature that preserves entropy per unit mass at each pressure level.
It writes the profile, diagnostics and plots. It does not use the UV solver.
Run the example
From the repository root, with the current kintera extension installed and NumPy, PyYAML, Torch and Matplotlib available:
python examples/example_jupiter_adiabat.py --output-dir examples/output/jupiter_adiabat
pytest -q tests/test_jupiter_tp_adiabat.py
Optional arguments are --yaml PATH, --levels N and --output-dir PATH.
The default card spans 6000 to 200 bar, starting at 2300 K, on 121 logarithmically
spaced pressure levels. This domain keeps the resulting temperatures inside
800–2500 K. It covers the two deep cloud transitions, not the visible weather
layers or the entire planetary interior. The boundary state and abundances
are illustrative inputs, not a fit to observations.
Reversible TP column. The clouds are retained and H2 evolves through the MgSiO3 reaction. Zero mass fractions are displayed at the plotting floor; no species floor is applied to the physical inventories.
Outputs and interpretation
profile.csv contains pressure, temperature, height above the lower boundary,
bulk density, specific entropy, iteration counts, saturation ratios and
conservation errors. x_NAME is the mole fraction of all gas and condensed
species combined; it is not a gas-only mixing ratio. q_NAME is the mass
fraction of the total parcel. summary.json records maximum residuals and
the first sampled cloudy level. profile.svg and profile.png contain
the three-panel plot.
Quantity |
Value |
|---|---|
First sampled MgSiO3 cloud level |
4519.16 bar; 2114.25 K |
First sampled KCl cloud level |
314.760 bar; 957.561 K |
Maximum specific-entropy error |
9.99e-6 J kg^-1 K^-1 |
Maximum relative element-inventory error |
2.85e-15 |
Maximum relative mass error |
2.23e-15 |
Maximum log saturation/complementarity residual |
9.73e-9 |
Cloud locations are grid samples at a 1e-12 cloud mole-fraction threshold, not a refined root of the cloud-base condition. The reported error values are numerical diagnostics for this approximation; they are not uncertainties on Jupiter’s cloud locations.
Complete profile CSVRecorded diagnostics JSON
Thermodynamic construction
At each new pressure, a bracketed scalar search calls TP equilibrium and solves
The density and entropy include the retained solids. Conserving specific entropy is necessary here: MgSiO3 formation changes the total number of molecules, so holding entropy per mole fixed would describe a different path. All temperature trials at one pressure restart from the same preceding parcel composition, preserving its elemental inventory. The first He species is an inert amount reference; H2 is a reacting gas, not a fixed reservoir.
For the ideal-gas card, concentrations are calculated directly as
This preserves trace SiH4 below the general TPX->V conversion’s default
1e-20 mol/m³ floor. Such trace concentrations still matter in the reaction
quotient. The same concentration scale applies to the suspended condensates.
Native TPV->S supplies the entropy, and relative_humidity(..., ngas)
checks the signed gas quotient including the five H2 product factors.
Altitude is a diagnostic obtained from hydrostatic balance with the specified constant gravity:
The script integrates this equation by the trapezoidal rule on the pressure grid. Temperature is found from entropy directly, rather than integrating an assumed dry lapse rate. The tests check 31/61-level agreement at shared pressure levels and recover the analytic constant-cp dry adiabat when all condensable inventories are removed.
Data choices and scope
The KCl curve is the registered kcl_lodders function: the bar-pressure fit
in Morley et al. (2012), Eq. 18, converted
to natural logarithms and Pa. The script checks that a KCl solid cloud never
appears above 1044 K. Evaluating the curve at cloud-free hotter levels is a
mathematical continuation, not a model of liquid KCl.
For MgSiO3 the example deliberately uses the conditional reconstruction documented in Magnesium silicate, from Visscher et al. (2010), Tables 2 and 3. It is supplied inline in the YAML; no unverified named formula is added to the production registry. The expression is
Its reconstruction assumes an H2-dominated mixture with X_H2 approximately 0.84. The explicit H2 inventory in this example is allowed to evolve. The full reaction quotient has net pressure exponent zero. The generic Antoine evaluator adds ln(100000), so the YAML subtracts 5 from its A coefficient to represent this dimensionless quotient correctly.
Heat capacities are an internally consistent idealized constant-cp closure, not measured species thermochemistry. For both curves L = ln(10)(A - B/T), the card enforces
The solid energy offsets and heat capacities are chosen accordingly, so the latent heat used in the entropy calculation agrees with the temperature derivative of the equilibrium curve. H2 rotational/vibrational variation, nonideal EOS effects at thousands of bars, and competing silicates such as forsterite are omitted. This example demonstrates the coupled numerical method; quantitative Jupiter predictions need a broader reaction network, verified caloric data, and an appropriate nonideal EOS.
YAML configuration
# Reduced, ideal-gas deep-Jupiter example, not a calibrated atmosphere.
# Reversible ascent: both condensates remain suspended in the parcel.
# Constant heat capacities below are an idealized closure, not measured data.
# They satisfy delta(cp/R)=0 for each constant-A-B/T equilibrium curve.
reference-state: {Tref: 0., Pref: 1.e5}
dynamics:
equation-of-state: {max-iter: 200, ftol: 1.e-8}
species:
- {name: He, phase: gas, composition: {He: 1}, cv_R: 1.5}
- {name: H2, phase: gas, composition: {H: 2}, cv_R: 2.5}
- {name: KCl, phase: gas, composition: {K: 1, Cl: 1}, cv_R: 2.5}
- {name: Mg, phase: gas, composition: {Mg: 1}, cv_R: 1.5}
- {name: SiH4, phase: gas, composition: {Si: 1, H: 4}, cv_R: 8.0}
- {name: H2O, phase: gas, composition: {H: 2, O: 1}, cv_R: 3.0}
# u0/R = -B*ln(10); with delta cp=0 this gives delta h/R=-B*ln(10).
- {name: KCl(s), phase: solid, composition: {K: 1, Cl: 1}, cv_R: 3.5, u0_R: -26208.02352845823}
- {name: MgSiO3(s), phase: solid, composition: {Mg: 1, Si: 1, O: 3}, cv_R: 6.0, u0_R: -114887.48321493792}
reactions:
# Morley et al. (2012), Eq. 18: log10(p_KCl/bar)=7.611-11382/T.
# https://arxiv.org/abs/1206.4313
- equation: KCl <=> KCl(s)
type: nucleation
rate-constant: {formula: kcl_lodders}
# CONDITIONAL reconstruction from Visscher et al. (2010), Tables 2 and 3:
# https://arxiv.org/html/1001.3639
# log10 Q = 5.86 - 5*log10(0.84) - 49895/T, Q=pMg*pSiH4*pH2O^3/pH2^5.
# This is NOT a directly published general reaction fit. See docs/reactions/mgsio3.
# The generic Antoine evaluator adds log10(1e5), hence A=6.238603569690592-5.
# Net pressure exponent is zero, so Q is dimensionless and no Pa/bar shift applies.
- equation: Mg + SiH4 + 3 H2O <=> MgSiO3(s) + 5 H2
type: nucleation
rate-constant: {formula: antoine, A: 1.238603569690592, B: 49895., C: 0.}
# Consumed by example_jupiter_adiabat.py; ignored by the native options reader.
column:
bottom_pressure_bar: 6000.
bottom_temperature_K: 2300.
top_pressure_bar: 200.
levels: 121
gravity_m_s2: 24.79 # Assumed constant gravity; ideal-gas hydrostatic column.
entropy_tolerance_J_kg_K: 1.e-5
temperature_bounds_K: [800., 2500.]
# Illustrative inventory with Mg/Si=1.2; not a fit to measured abundances.
mole_fractions:
He: 0.14
H2O: 0.001
Mg: 0.00006
SiH4: 0.00005
KCl: 0.0000002
H2: balance
Python implementation
#!/usr/bin/env python3
"""Reversible 1D deep-Jupiter adiabat with KCl and MgSiO3 condensation.
Each temperature trial calls ThermoX.forward (TP equilibrium). A bracketed
outer solve conserves specific entropy, including retained condensates.
Run: python examples/example_jupiter_adiabat.py --output-dir /tmp/jupiter
The MgSiO3 curve and constant heat capacities are illustrative approximations.
"""
import argparse
import csv
import json
from pathlib import Path
import numpy as np
import torch
import yaml
from kintera import ThermoOptions, ThermoX, relative_humidity, constants
DEFAULT_CARD = Path(__file__).with_name('jupiter_kcl_mgsio3.yaml')
def run_column(card=DEFAULT_CARD, levels=None):
"""Return an equilibrium isentrope and independent conservation diagnostics."""
config = yaml.safe_load(Path(card).read_text())
settings = config['column']
levels = levels if levels is not None else settings['levels']
if levels < 2:
raise ValueError('At least two pressure levels are required')
pb = float(settings['bottom_pressure_bar']) * 1e5
pt = float(settings['top_pressure_bar']) * 1e5
if not 0 < pt < pb:
raise ValueError('Require 0 < top pressure < bottom pressure')
tmin, tmax = map(float, settings['temperature_bounds_K'])
tbase = float(settings['bottom_temperature_K'])
if not 0 < tmin <= tbase <= tmax:
raise ValueError('Bottom temperature must be inside the positive bounds')
gravity = float(settings['gravity_m_s2'])
tolerance = float(settings['entropy_tolerance_J_kg_K'])
if gravity <= 0 or tolerance <= 0:
raise ValueError('Gravity and entropy tolerance must be positive')
thermo = ThermoX(ThermoOptions.from_yaml(str(card)))
names = thermo.options.species()
ngas = len(thermo.options.vapor_ids())
he = names.index('He') # Inert amount reference; H2 participates in reaction.
mole = settings['mole_fractions']
unknown = set(mole) - set(names)
if unknown:
raise ValueError(f'Unknown initial species: {unknown}')
x0 = np.array([float(mole.get(n, 0.)) if mole.get(n) != 'balance' else 0.
for n in names])
balance = [i for i, name in enumerate(names) if mole.get(name) == 'balance']
if len(balance) != 1:
raise ValueError('Exactly one species must have mole fraction balance')
x0[balance[0]] = 1. - x0.sum()
if np.any(x0 < 0) or x0[he] <= 0:
raise ValueError('Initial fractions must be nonnegative with positive He')
seed = torch.tensor(x0[None, :], dtype=torch.float64)
mu = thermo.buffer('mu').numpy()
compositions = {s['name']: s['composition'] for s in config['species']}
elements = sorted({e for c in compositions.values() for e in c})
atoms = np.array([[compositions[n].get(e, 0.) for n in names] for e in elements])
inventory = atoms @ x0
mass0 = x0 @ mu
stoich = thermo.buffer('stoich')
pressures = np.geomspace(pb, pt, levels)
rows, previous = [], seed
def equilibrate(temperature, pressure, start):
# Every trial starts from the same inventory, never the preceding
# temperature trial. That avoids accumulating root-search history.
x = start.clone()
t = torch.tensor([temperature], dtype=torch.float64)
p = torch.tensor([pressure], dtype=torch.float64)
diag = torch.zeros((1, 1), dtype=torch.float64)
thermo.forward(t, p, x, False, diag)
if diag.item() < 0 or not torch.isfinite(x).all() or (x < 0).any():
raise RuntimeError(f'TP failed at {temperature:g} K, {pressure/1e5:g} bar: {diag.item()}')
# This card is ideal gas. Preserve depleted SiH4 below the general
# TPX->V conversion's 1e-20 mol/m3 floor when checking equilibrium.
conc = x * p[:, None] / (constants.Rgas*t[:, None]*x[:, :ngas].sum(-1, keepdim=True))
rho = (conc * thermo.buffer('mu')).sum(-1)
entropy = thermo.compute('TPV->S', [t, p, conc]) / rho
return x, conc, float(rho.item()), float(entropy.item()), int(diag.item())
with torch.no_grad():
base = equilibrate(tbase, pb, seed)
target = base[3]
for level, pressure in enumerate(pressures):
if level == 0:
temperature, state, nouter = tbase, base, 0
else:
# An ascent cools the parcel. Use the previous equilibrium as
# the inventory-preserving starting composition for all trials.
lo, hi = max(tmin, rows[-1]['temperature_K']*.65), rows[-1]['temperature_K']
low = equilibrate(lo, pressure, previous)
high = equilibrate(hi, pressure, previous)
if not low[3] <= target <= high[3]:
raise RuntimeError('Entropy root is outside temperature bounds; '
'raise the top pressure or revise the model domain')
for nouter in range(1, 65):
temperature = .5*(lo+hi)
state = equilibrate(temperature, pressure, previous)
if abs(state[3]-target) <= tolerance:
break
if state[3] < target:
lo = temperature
else:
hi = temperature
else:
raise RuntimeError('Specific-entropy root failed to converge')
x, conc, rho, entropy, diag = state
amounts = x[0].numpy() * x0[he] / x[0, he].item()
atom_scale = np.where(inventory > 0, inventory, 1.)
atom_error = float(np.max(np.abs(atoms @ amounts-inventory)/atom_scale))
mass_error = float(abs(amounts @ mu / mass0-1))
rh = relative_humidity(torch.tensor([temperature], dtype=torch.float64),
conc, stoich, thermo.options.nucleation(), ngas)[0].numpy()
cloud = x[0, ngas:].numpy()
with np.errstate(divide='ignore'):
violation = np.where(cloud > 1e-15, abs(np.log(rh)), np.maximum(0, np.log(rh)))
if not np.isfinite(rh).all() or np.max(violation) > 2e-7:
raise RuntimeError(f'Phase equilibrium check failed at level {level}: T={temperature}, diag={diag}, RH={rh}, x={x}')
row = dict(pressure_bar=pressure/1e5, temperature_K=temperature,
density_kg_m3=rho, entropy_J_kg_K=entropy,
entropy_error_J_kg_K=entropy-target,
element_relative_error=atom_error, mass_relative_error=mass_error,
tp_iterations=diag, entropy_iterations=nouter,
rh_KCl=float(rh[0]), rh_MgSiO3=float(rh[1]),
phase_residual=float(np.max(violation)))
if level == 0:
row['height_km'] = 0.
else:
# dP/dz=-rho*g, with bulk density including suspended clouds.
h0 = pressures[level-1]/rows[-1]['density_kg_m3']/gravity
h1 = pressure/rho/gravity
row['height_km'] = rows[-1]['height_km']-.5*(h0+h1)*np.log(pressure/pressures[level-1])/1000
for i, name in enumerate(names):
row['x_'+name] = float(x[0, i]) # Fraction of all gas+cloud moles.
row['q_'+name] = float(x[0, i]*mu[i]/(x[0].numpy()@mu))
rows.append(row); previous = x
# Only use the solid KCl branch where a cloud is actually present.
if any(r['x_KCl(s)'] > 1e-15 and r['temperature_K'] >= 1044 for r in rows):
raise RuntimeError('Solid KCl cloud encountered above its melting point')
if max(r['element_relative_error'] for r in rows) > 1e-7:
raise RuntimeError('Element conservation check failed')
summary = dict(levels=levels, species=names, closure='reversible; condensates retained',
MgSiO3_fit='conditional Visscher table reconstruction; illustrative',
max_entropy_error_J_kg_K=max(abs(r['entropy_error_J_kg_K']) for r in rows),
max_element_relative_error=max(r['element_relative_error'] for r in rows),
max_mass_relative_error=max(r['mass_relative_error'] for r in rows),
max_phase_residual=max(r['phase_residual'] for r in rows))
for name in ('KCl(s)', 'MgSiO3(s)'):
cloudy = [r for r in rows if r['x_'+name] > 1e-12]
summary[name+'_first_cloudy_level'] = ({k: cloudy[0][k] for k in ('pressure_bar','temperature_K','height_km')}
if cloudy else None)
return rows, summary
def write_outputs(rows, summary, output):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
output = Path(output); output.mkdir(parents=True, exist_ok=True)
with (output/'profile.csv').open('w') as file:
writer = csv.DictWriter(file, fieldnames=list(rows[0]), lineterminator='\n')
writer.writeheader(); writer.writerows(rows)
(output/'summary.json').write_text(json.dumps(summary, indent=2)+'\n')
p = np.array([r['pressure_bar'] for r in rows])
fig, axes = plt.subplots(1, 3, figsize=(12, 5), sharey=True)
axes[0].plot([r['temperature_K'] for r in rows], p)
axes[0].set_xlabel('Temperature [K]'); axes[0].set_ylabel('Pressure [bar]')
for name in ('KCl', 'KCl(s)', 'Mg', 'SiH4', 'MgSiO3(s)'):
q = np.array([r['q_'+name] for r in rows])
axes[1].plot(np.maximum(q, 1e-18), p, label=name)
axes[1].set_xscale('log'); axes[1].set_xlim(1e-12, 1e-2)
axes[1].set_xlabel('Mass fraction'); axes[1].legend(fontsize=8)
for name in ('KCl', 'MgSiO3'):
axes[2].plot([r['rh_'+name] for r in rows], p, label=name)
axes[2].axvline(1, color='0.5', linestyle=':')
axes[2].set_xscale('log'); axes[2].set_xlabel('Saturation ratio Q / Qsat')
axes[2].legend()
for ax in axes:
ax.set_yscale('log'); ax.grid(alpha=.25)
axes[0].invert_yaxis()
fig.suptitle('Illustrative deep-Jupiter reversible TP adiabat')
fig.tight_layout()
plt.rcParams['svg.hashsalt'] = 'jupiter-tp-adiabat'
fig.savefig(output/'profile.svg', metadata={'Date': None})
svg = output/'profile.svg'
svg.write_text('\n'.join(line.rstrip() for line in svg.read_text().splitlines())+'\n')
fig.savefig(output/'profile.png', dpi=170)
plt.close(fig)
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--yaml', type=Path, default=DEFAULT_CARD)
parser.add_argument('--levels', type=int)
parser.add_argument('--output-dir', type=Path, default=Path('examples/output/jupiter_adiabat'))
args = parser.parse_args()
rows, summary = run_column(args.yaml, args.levels)
write_outputs(rows, summary, args.output_dir)
print(json.dumps(summary, indent=2))
print(f'Wrote profile.csv, summary.json, profile.svg and profile.png in {args.output_dir}')
if __name__ == '__main__':
main()