Numerical Methods

This page translates the numerical-methods part of the TARSA manuscript into documentation form. It summarizes how the governing transport equation is discretized and advanced in time.

Grid, indexing, and state vector

The transport equation is solved on a structured finite-volume grid with $N_x \times N_y \times N_z$ cells. Each cell $(i,j,k)$ has dimensions $\Delta x_i$, $\Delta y_j$, and $\Delta z_k$. Horizontal cell sizes depend on the chosen projection, while vertical cell thicknesses are derived from the ERA5-based height coordinate.

For efficient implementation, the three-dimensional grid is mapped to a one-dimensional state vector using Julia column-major storage for arrays ordered (Nx, Ny, Nz): i is fastest, then j, and vertical level k is slowest. The index mapping is

\[p = (k - 1)N_xN_y + (j - 1)N_x + i,\]

with $p = 1, \dots, N_xN_yN_z$. The state vector at time level $n$ is therefore

\[\mathbf{X}^n = \left(C_1^n, C_2^n, \dots, C_N^n\right)^{\mathrm T}.\]

Finite-volume discretization

To keep the scheme mass consistent, TARSA discretizes the two conservative balances from the governing equations — the air-mass budget and the tracer-mass budget — over each control volume $V_{i,j,k}$. Integrating both in flux form gives

\[V_{i,j,k}\frac{d\rho_{i,j,k}}{dt} + \int_{\partial V_{i,j,k}} (\rho\mathbf{u})\cdot\mathbf{n}\,dS = 0,\]

\[V_{i,j,k}\frac{d(\rho C)_{i,j,k}}{dt} + \int_{\partial V_{i,j,k}} (\rho\mathbf{u}C)\cdot\mathbf{n}\,dS = \int_{\partial V_{i,j,k}} \rho K_z \frac{\partial C}{\partial z}\,\mathbf{e}_z\cdot\mathbf{n}\,dS + V_{i,j,k}\,\rho_{i,j,k}E_{i,j,k},\]

where $\mathbf{n}$ is the outward unit normal. In the discrete system the surface integrals are replaced by flux sums over the six cell faces. Each face carries an air-mass flux $G$ and a tracer-mass flux $G_\varphi$, built from the same reconstructed face winds and densities, so the two budgets stay tied together. The mixing ratio is then recovered as $C = (\rho C)/\rho$. This discrete pairing is what guarantees a uniform $C$ is preserved exactly (the constancy invariant); see the advection section below.

Advection schemes

TARSA provides two finite-volume advection options. The first, denoted $\mathcal{A}_{\mathrm{UP}}$, is a robust first-order upwind scheme; it is monotone by construction and is used for the monolithic fully implicit CAMS replay configuration. The second, denoted $\mathcal{A}_{\mathrm{rKM}}$, is the regularized Koren MUSCL scheme — a higher-order TVD-limited piecewise-linear MUSCL reconstruction with the cubic-Hermite-regularized Koren limiter — and is the default for the sharp-plume experiments. In the code this is the mass-consistent DST-Koren family (scheme id :dst_koren_rk3), which advects the density-weighted pair $(\rho, \rho C)$; the first-order upwind scheme is a single-field volume-flux update used by the implicit replay solver.

First-order upwind ($\mathcal{A}_{\mathrm{UP}}$)

The face-normal advective flux is

\[F_{\mathrm{face}} = (\mathbf{u}_{\mathrm{face}}\cdot\mathbf{n}_{\mathrm{face}})S_{\mathrm{face}},\]

where $S_{\mathrm{face}}$ is the face area. At a face shared by adjacent cells, the transported value is taken from the upwind cell according to the sign of the face flux.

This option is the most diffusive of the available schemes, but it is robust and is used by the fully implicit replay solver. It is a single-field (volume-flux) update and is not constancy-preserving; mass consistency in TARSA is provided by the explicit DST-Koren family described next, which is the default.

Regularized Koren MUSCL mass-consistent explicit advection ($\mathcal{A}_{\mathrm{rKM}}$)

For sharp-plume calculations, TARSA uses a higher-order TVD-limited explicit scheme based on piecewise-linear MUSCL reconstruction with a Koren-family limiter, applied in mass-consistent form: instead of advecting the mixing ratio $C$ alone, the scheme advances the density-weighted pair $(\rho, \rho C)$ with a single set of face mass fluxes and recovers $C = \rho C/\rho$ at the end of the step. Because the limiter depends on local solution gradients, $\mathcal{A}_{\mathrm{rKM}}$ is nonlinear at the discrete level (its one-step Jacobian is trajectory-dependent), unlike the affine $\mathcal{A}_{\mathrm{UP}}$ update.

The Koren limiter itself is regularized with a $C^1$ cubic-Hermite smoothing near the classical piecewise breakpoints. This removes limiter kinks while preserving the TVD admissible region, which is useful for gradient-based applications. At each face, the same smoothed limiter is applied independently to $C$ and to $\rho$, giving Koren-limited left/right interface states $C_L, C_R$ and $\rho_L, \rho_R$. These are combined through the smooth flux split

\[u_{\delta}^{\pm} = \frac{1}{2}\left(u \pm \sqrt{u^2 + \delta^2}\right),\]

with a small smoothing parameter $\delta > 0$, to form the two face fluxes

\[G_{\mathrm{face}} = S_{\mathrm{face}}\left(u_{\delta}^{+}\rho_L + u_{\delta}^{-}\rho_R\right), \qquad G_{\varphi,\mathrm{face}} = S_{\mathrm{face}}\left(u_{\delta}^{+}\rho_L C_L + u_{\delta}^{-}\rho_R C_R\right),\]

where $G$ is the air-mass flux and $G_\varphi$ the tracer-mass flux. The cell tendencies are the discrete divergences

\[\frac{d\rho_{i,j,k}}{dt} = -\frac{1}{V_{i,j,k}}\sum_{\mathrm{faces}} G_{\mathrm{face}}, \qquad \frac{d(\rho C)_{i,j,k}}{dt} = -\frac{1}{V_{i,j,k}}\sum_{\mathrm{faces}} G_{\varphi,\mathrm{face}},\]

after which $C = \rho C/\rho$ is recovered.

Constancy invariant. When $C \equiv c$ everywhere, the interface states satisfy $C_L = C_R = c$, so $G_{\varphi,\mathrm{face}} = c\,G_{\mathrm{face}}$ on every face. The tracer-mass update is then exactly $c$ times the air-mass update, $\rho C^\star = c\,\rho^\star$, and the recovered $C^\star = c$ — to floating-point roundoff, regardless of how divergent the discrete wind is. The density $\rho$ is a prescribed parameter, frozen as an input each step; it is "transported" only so that its discrete divergence matches the one applied to $\rho C$, which is precisely what makes the cancellation exact.

Boundary conventions. The density $\rho$ uses zero-gradient (ghost = nearest interior cell) conditions at the lateral walls and at both the model bottom and the open model top, so a uniform field passes through the top boundary unchanged. The mixing ratio $C$ retains its own boundary treatment (Dirichlet inflow or zero-gradient) independently of $\rho$.

Because this advection step is explicit, the time step must satisfy a CFL condition:

\[\Delta t \le \frac{\mathrm{CFL}}{ \max\limits_{i,j,k} \left( \frac{|u|}{\Delta x_i} + \frac{|v|}{\Delta y_j} + \frac{|w|}{\Delta z_k} \right) }.\]

Representative sharp-plume calculations use $\mathrm{CFL} \approx 0.5$.

Vertical diffusion

Vertical diffusive fluxes are discretized with a second-order central-difference formulation. Vertical gradients are approximated at cell interfaces, and both $K_z$ and $\rho$ are interpolated to those interfaces. Because the vertical grid is generally non-uniform, TARSA uses distance-weighted linear interpolation at faces, which reduces to the arithmetic mean on a uniform grid.

The resulting implicit vertical-diffusion operator couples only neighboring levels within each column.

Time integration and solver configurations

Different experiments use different solver combinations, but the building blocks are the same:

• explicit SSP-RK3 for advection
• implicit column solves for vertical diffusion
• second-order Strang splitting to combine advection and diffusion
• an optional monolithic fully implicit backward-Euler configuration for replay problems with strongly constrained boundaries

SSP-RK3 update

In the mass-consistent scheme the integrated state is the conservative pair $\mathbf{W} = (\rho, \rho C)$, and $L(\mathbf{W})$ denotes the Koren mass-flux divergence operator that returns both tendencies $(-\nabla\!\cdot G, -\nabla\!\cdot G_\varphi)$ from the advection section. The strong-stability-preserving third-order Runge-Kutta update is

\[\mathbf{W}^{(1)} = \mathbf{W}^{n} + \Delta t\,L\bigl(\mathbf{W}^{n}\bigr),\]

\[\mathbf{W}^{(2)} = \frac{3}{4}\,\mathbf{W}^{n} + \frac{1}{4}\left( \mathbf{W}^{(1)} + \Delta t\,L\bigl(\mathbf{W}^{(1)}\bigr) \right),\]

\[\mathbf{W}^{n+1} = \frac{1}{3}\,\mathbf{W}^{n} + \frac{2}{3}\left( \mathbf{W}^{(2)} + \Delta t\,L\bigl(\mathbf{W}^{(2)}\bigr) \right),\]

after which $C^{n+1} = (\rho C)^{n+1}/\rho^{n+1}$ is recovered. Because every RK3 stage applies the same convex combination to both components of $\mathbf{W}$, the constancy invariant survives the multi-stage update: $C \equiv c$ implies $(\rho C)^{(s)} = c\,\rho^{(s)}$ at every stage. The RK3 right-hand side evaluates all three directions together, so for a divergence-free flow $\rho$ is left unchanged ($-\nabla\!\cdot(\rho\mathbf{u}) = 0$).

Implicit vertical diffusion

The vertical diffusion step leads to an independent tridiagonal solve in each $(i,j)$ column:

\[\mathbf{T}_{i,j}\,\mathbf{X}_{i,j}^{n+1} = \mathbf{r}_{i,j}.\]

Each system is solved with the Thomas algorithm in $\mathcal{O}(N_z)$ operations per column.

Strang splitting

In the default sharp-plume configuration, advection and vertical diffusion are combined with second-order Strang splitting:

\[\mathbf{X}^{n+1} = \Phi_{\mathrm{diff}}(\Delta t/2)\, \Phi_{\mathrm{adv}}(\Delta t)\, \Phi_{\mathrm{diff}}(\Delta t/2)\, \mathbf{X}^{n} + \mathbf{q},\]

where $\mathbf{q}$ collects source, sink, and boundary contributions over the split step.

The DST-Koren advection sub-step itself offers a directional (dimensional) split variant as an alternative to SSP-RK3 (the :dst_koren_dimsplit stepper, distinct from the outer Strang process split above): the conservative pair is swept $x(\tfrac{\Delta t}{2}) \to y(\tfrac{\Delta t}{2}) \to z(\Delta t) \to y(\tfrac{\Delta t}{2}) \to x(\tfrac{\Delta t}{2})$, with each one-dimensional sweep applying the same density-weighted Koren fluxes to both $\rho$ and $\rho C$. Carrying the pair through every sweep keeps the constancy invariant exact. The directional split does introduce an $O(\Delta t^2)$ drift in the internal working $\rho$ relative to the prescribed input density; the recovered $C = \rho C/\rho$ is unaffected (numerator and denominator carry the same drift), but an external mass diagnostic that pairs $\rho C$ with the fixed input $\rho$ is exact only for the RK3 stepper, which evaluates all directions simultaneously.

Fully implicit replay option

For replay-style experiments with strongly prescribed, time-varying boundaries, TARSA also provides a monolithic fully implicit backward-Euler configuration with first-order upwind advection and implicit vertical diffusion assembled into a global sparse linear system per time step. This configuration is more diffusive than the sharp-plume setup and uses the single-field volume-flux advection, so it is not constancy-preserving; it is retained because it is numerically robust for constrained replay problems such as the CAMS consistency experiment, where the prescribed boundaries dominate the interior solution.

Typical experiment configurations

The manuscript uses the following solver combinations. All regularized Koren MUSCL ($\mathcal{A}_{\mathrm{rKM}}$) rows use the mass-consistent $(\rho, \rho C)$ formulation (the default, mass_consistent=true, scheme id :dst_koren_rk3); only the first-order-upwind ($\mathcal{A}_{\mathrm{UP}}$) CAMS replay row is non-mass-consistent.