Basic equations and numerics

VULCAN solves the one-dimensional Eulerian continuity equation (Section 2.1 of Tsai et al. 2021 ¹) for the number density of each species. This page covers the transport flux, its spatial discretization, and the time-integration scheme.

Continuity equation

For each species \(\text{i}\),

with \(n_{\text{i}}\) the number density (cm\(^{-3}\)), \(P_{\text{i}}\) and \(L_{\text{i}}\) the chemical production and loss rates (cm\(^{-3}\) s\(^{-1}\)), and \(\Phi_{\text{i}}\) the vertical transport flux. \(P_{\text{i}} - L_{\text{i}}\) is evaluated by chemdf in the generated module chem_funs.py; see chemical networks.

Transport flux

Assuming hydrostatic balance, the flux combines advection, eddy diffusion, and molecular plus thermal diffusion (following Chamberlain & Hunten 1987 ³):

where \(v\) is the vertical wind, \(K_{zz}\) the eddy diffusion coefficient, \(D_{\text{i}}\) the molecular diffusion coefficient, \(X_{\text{i}} = n_{\text{i}}/n_\mathrm{tot}\) the mixing ratio, \(H_{\text{i}} = k_B T/(m_{\text{i}} g)\) the species scale height, and \(\alpha_T\) the thermal diffusion factor. The four contributions are, in order: advection along the wind, eddy diffusion smearing out compositional gradients, molecular diffusion driving each species toward its own diffusive equilibrium, and thermal diffusion (a secondary effect except for light species in the thermosphere).

Spatial discretization

The diffusive terms use a second-order central difference and the advective term a first-order upwind scheme (Brasseur & Jacob 2017 ⁴). On the staggered grid the flux divergence in layer \(\text{j}\) is

Substituting Eq. (2) reduces Eq. (1) to a system of ODEs of the form \(A_{\text{j}} y_{\text{j}} + B_{\text{j}+1} y_{\text{j}+1} + C_{\text{j}-1} y_{\text{j}-1}\) per layer.

Implementation in VULCAN

The transport operator is assembled in op.ODESolver, with four variants selected at runtime in Ros2.solver according to the config flags:

The upwind advection appears as the (vz>0)/(vz<0) masks. Interface quantities (\(T_{\text{i}}\), \(H_{p,\text{i}}\), \(\Delta z_{\text{i}}\)) are precomputed in build_atm.Atm.f_mu_dz, and the molecular diffusion coefficients \(D_{zz}\) and thermal diffusion factors \(\alpha_T\) in build_atm.Atm.mol_diff. The latter scales a reference binary diffusion coefficient (per background gas atm_base: H\(_2\), N\(_2\), O\(_2\), or CO\(_2\), after Marrero & Mason 1972 ⁶ / Banks & Kockarts 1973 ⁷) by molecular mass, corresponding to Appendix A of the paper ¹.

Time integration

The stiff ODE system is integrated with a second-order Rosenbrock method ⁵ in op.Ros2.solver. For \(\mathrm{d}n/\mathrm{d}t = f(n)\), explicit (forward-Euler) differencing is cheap but unstable at the large steps stiff chemistry needs, while fully implicit (backward-Euler) differencing is stable but requires an expensive Newton–Raphson iteration for \(n_{\text{k}+1}\). The Rosenbrock method is linearly implicit: it uses the Jacobian \(J \equiv \partial f/\partial n\) directly, avoiding iteration. The two-stage, second-order scheme is

with \(\gamma = 1 + 1/\sqrt{2}\). Verwer et al. (1999) ⁵ recommend this method for its stability over large stepsizes, which suits the needs of chemical kinetics. Both stages share the same left-hand-side matrix

so it is factored only once per step. The chemical part of \(J\) is analytic (generated symbolically by make_chem_funs.make_jac/make_neg_jac into chem_funs.symjac/neg_symjac); the transport part is added in the lhs_jac_* methods. The combined block-tridiagonal matrix is stored in banded form (store_bandM) and solved with scipy.linalg.solve_banded.

Step-size control and convergence

The step size adapts to the truncation error \(\mathcal{E} = |n_{\text{k}+1} - n^*_{\text{k}+1}|\), where \(n^*_{\text{k}+1} = n_{\text{k}} + \Delta t\, g_1\) is the embedded first-order estimate. In Ros2.step_size,

where rtol is the desired relative tolerance and 0.9 is a safety factor. The step-change factor is bounded to \([r_{\min}, r_{\max}] = [0.5,\ 2]\) — the config parameters dt_var_min and dt_var_max — and the step itself to \([10^{-8},\ 10^{18}]\) s (dt_min, dt_max). A step with \(\mathcal{E} > \mathrm{rtol}\) is rejected and retried with a smaller stepsize.

Steady state is declared in op.Integration from the long-term relative variation, following the criteria of Tsai et al. (2017) ², evaluated over a look-back window from \(f\tau\) to the current integration time \(\tau\) (with \(f\) set by st_factor):

where \(\text{k}'\) is the timestep at \(f\tau\). Species with mixing ratios below mtol_conv or number densities below atol are excluded from the maximum so that negligible trace species do not trigger it spuriously. Convergence requires both \(\Delta\hat{n}\) and its rate \(\Delta\hat{n}/\Delta t\) to fall below either the primary thresholds (yconv_cri, slope_cri) or the relaxed thresholds (yconv_min, slope_min), where the relaxed slope adapts to the local vertical-diffusion timescale, \(\sim \min_{\text{j}} K_{zz,\text{j}}/(0.1\,H_{p,\text{j}})^2\). With photochemistry enabled, both conditions additionally require the actinic-flux change to fall below flux_cri.