TChem-atm (v2.0.0): scalable performance-portable multiphase atmospheric chemistry

We present TChem-atm, a performance-portable approach that enables efficient simulation of chemically detailed and multiphase atmospheric chemistry on modern heterogeneous computing architectures. Unlike previous efforts that rely on architecture-specific code or focus exclusively on gas-phase chemistry, TChem-atm supports fully coupled gas–aerosol systems with execution across CPUs, NVIDIA GPUs, and AMD GPUs through the Kokkos programming model. It integrates the flexible multiphase capabilities of the Community Atmospheric Model Chemistry Package (CAMP) with the high-performance kinetic routines of TChem, and includes automatic Jacobian construction with support for a range of stiff ODE solvers. In a proof-of-concept integration with the particle-resolved model PartMC, TChem-atm reproduces the existing PartMC–CAMP implementation within solver tolerances and delivers substantial GPU speedups, especially for large particle populations. Performance benchmarks reveal substantial speedups on GPU platforms, particularly for large particle populations, with consistent results across hardware backends. TChem-atm enables performance-portable execution across CPUs and GPUs, though optimal efficiency may require modest architecture-specific tuning (e.g., team and vector sizes), with up to a twofold improvement on the NVIDIA H100. It directly supports sectional and particle-resolved host models, while modal aerosol schemes require minor adaptation to provide particle-scale quantities such as representative diameters. By enabling chemically detailed, multiphase simulations with performance portability and host-model flexibility, TChem-atm facilitates the incorporation of advanced chemistry into atmospheric models.

1 Introduction

Simulating the chemical evolution of the atmosphere is central to predicting air quality, climate forcing, and ecosystem impacts. Atmospheric models must capture complex interactions between gas-phase reactions, aerosol microphysics, and multiphase chemistry – often at fine spatial and temporal resolutions and over extended simulation domains. As chemical mechanisms become more detailed and models increase in resolution, the computational cost of solving the associated systems of equations grows rapidly. To address this challenge, flexible and scalable methods are needed that can integrate complex chemistry into models while taking full advantage of modern computing architectures.

Simulating detailed gas-phase mechanisms, multiphase processes, and aerosol microphysics at such fine resolutions typically requires solving thousands of coupled, numerically stiff ordinary differential equations (ODEs) at each grid point. This makes chemical solvers one of the most computationally intensive components of atmospheric models. Graphics Processing Units (GPUs) offer a massively parallel architecture that can accelerate these calculations significantly, enabling more detailed and efficient simulations.

This acceleration is critical across a wide range of applications. In Earth system models such as E3SM (Energy Exascale Earth System Model), interactive chemistry and aerosols are essential for representing radiative forcing and cloud feedbacks over long timescales (Golaz et al., 2019; Rasch et al., 2019). In global chemical transport models like GEOS-Chem, GPU acceleration enables higher spatial resolution and real-time data assimilation (Bey et al., 2001; Eastham et al., 2018). Operational air quality models such as CMAQ must meet tight runtime constraints while still representing complex chemical and physical processes (Byun and Schere, 2006; Appel et al., 2021). Across these domains, reducing the computational burden of atmospheric chemistry is essential for improving both fidelity and performance.

Recent work has demonstrated the potential of GPU acceleration for atmospheric models. For example, Ruiz et al. (2024) achieved a 35× speed-up by distributing the CAMP solver load across GPU threads, and Alvanos and Christoudias (2017) reported kernel-level speed-ups of up to 20.4× for CUDA-based kinetic integration in EMAC. Cao et al. (2023) showed more than 1000× speed-ups for GPU-accelerated advection in CAMx using hybrid GPU–MPI strategies, while Quevedo et al. (2025) demonstrated that GPU-enabled CMAQ simulations could cut runtimes nearly in half. Similarly, Sun et al. (2018) highlighted the importance of memory layout and GPU-aware optimization in achieving efficient Rosenbrock solver performance in CAM4-Chem.

While these studies underscore the promise of GPU-accelerated atmospheric chemistry, most rely on architecture-specific APIs (e.g., CUDA) and hand-optimized kernels tightly coupled to the host model. This approach limits code portability, increases integration complexity, and hinders broader adoption of GPU acceleration across diverse computing environments.

To address these limitations, performance-portable approaches such as TChem have emerged, leveraging libraries like Kokkos to abstract architecture-specific details (Kim et al., 2023). TChem provides a scalable and efficient backend for evaluating gas-phase and surface reaction kinetics, thermodynamic properties, and Jacobians via automatic differentiation. However, TChem was originally developed for gas-phase chemistry and did not natively support multiphase processes, limiting its application in comprehensive atmospheric chemistry models.

In parallel, we developed the Chemistry Across Multiple Phases (CAMP) library as a runtime-configurable tool for simulating gas- and aerosol-phase chemistry in models with varying aerosol representations (Dawson et al., 2022). CAMP decouples chemical logic from the host model's aerosol structure and supports modular chemical mechanism configuration and solver abstraction. However, while CAMP is flexible in mechanism design, its reliance on host-side solvers limited its scalability and portability to new hardware architectures.

Recognizing their complementary capabilities, the present work combines CAMP's support for multiphase chemistry with TChem's performance-portable computational backend. We refer to this extended version as TChem-atm, which is maintained as a standalone code with its own repository. TChem-atm builds on the software infrastructure of TChem (v2.0), which was originally developed for gas-phase kinetics in combustion applications. TChem provides the foundational components that TChem-atm inherits, namely the kinetic-model data structures, automatic and numerical Jacobian construction, and the batched-evaluation framework implemented with Kokkos. TChem-atm extends this base by introducing atmospheric-chemistry capabilities not present in TChem, including (i) support for aerosol- and multiphase reaction mechanisms, (ii) integration with CAMP's gas–aerosol process abstractions, (iii) an aerosol-model constant-data object and associated particle/section/mode interfaces, and (iv) mechanism parsing for heterogeneous and condensed-phase chemistry. These additions transform TChem's gas-phase engine into a performance-portable multiphase chemistry library suitable for atmospheric models across CPUs and GPUs. TChem-atm further provides interoperability with multiple ODE solvers, including TINES (Kim and Diaz-Ibarra, 2021) and SUNDIALS (Balos et al., 2021; Gardner et al., 2022; Hindmarsh et al., 2005, 2025).

TChem-atm is compatible with a wide range of aerosol representations, but the degree of direct interoperability differs across host models. Sectional and particle-resolved frameworks map naturally onto TChem-atm's computational-particle abstraction. Modal aerosol schemes, however, typically do not expose particle-scale quantities such as representative diameters or per-species mass vectors, and therefore require modest adaptation to supply these inputs. Once provided, TChem-atm can evaluate multiphase chemistry within modal frameworks in the same manner as for other aerosol representations. While both CAMP and TChem provide flexibility in chemical mechanism configuration, CAMP offers a structure tailored for gas–aerosol interactions, whereas TChem enables efficient execution on heterogeneous architectures via the Kokkos programming model.

This paper presents the first application of this integrated library to atmospheric modeling, with a focus on enabling efficient aerosol–chemistry interactions in high-resolution simulations. The present study is intentionally scoped as a proof-of-concept implementation that validates correctness and establishes performance characteristics of the new multiphase chemistry framework. Because the integration of CAMP chemistry with TChem's performance-portable backend is novel, our immediate priority is demonstrating functional equivalence to existing implementations and quantifying solver and hardware scaling. Application within regional- or global-scale models is a natural next step, but lies beyond the scope of this initial paper. The ability to interchange solvers across platforms without reimplementing solver logic streamlines deployment and optimization. By supporting scalable chemistry calculations on diverse architectures, this work represents a concrete step toward the generalized aerosol/chemistry interface advocated by Hodzic et al. (2023).

Although TChem-atm achieves performance portability across CPUs and GPUs without architecture-specific code, performance is further improved when tunable parameters (e.g., Kokkos team and vector sizes) are optimized for the target hardware. As demonstrated in Sect. 3.2, architecture-aware tuning yields up to a factor-of-two speed improvement on the NVIDIA H100 relative to the Kokkos default configuration. We note that while our current implementation succeeds as a first step toward performance portability, individual applications could be tuned for even greater performance. Such tunings would depend on the chemistry mechanisms present in the application, the choice of solver, and the specific computational architecture. Here, we introduce the capability to flexibly define atmospheric chemistry mechanisms for a variety of applications in a common computational environment, capable of running on a variety of advanced computing architectures; optimized tunings for specific applications, solvers, and architectures is a subject marked for future work.

2 Methods

TChem-atm is a performance-portable chemistry library designed to support multiphase atmospheric chemistry across a range of computing architectures. Figure 1 presents a high-level overview of the TChem-atm library and outlines the structure of this section.

We begin with an overview of the software architecture and the relationship between TChem, CAMP, and TChem-atm (Sect. 2.1), since these design choices define the core abstractions used throughout the implementation. The workflow then begins with user-defined YAML input files that specify gas-phase chemical mechanisms and aerosol processes (Sect. 2.2). These inputs are parsed into structured data objects for the gas and aerosol phases (Sect. 2.3), which are synchronized to device memory via dual Kokkos views. TChem-atm then constructs source terms for gas-phase chemistry (Sect. 2.5) and for aerosol-phase processes using a flexible abstraction that supports sectional and particle-resolved representations (Sect. 2.7 and 2.6). Performance portability is achieved through the Kokkos programming model (Sect. 2.8), which enables execution on CPUs as well as NVIDIA and AMD GPUs. If needed, Jacobian matrices are computed (Sect. 2.9) and passed to one of several supported ODE solvers (Sect. 2.10). Our testing approach is described in Sect. 2.11. Finally, we describe the integration of TChem-atm with PartMC, which enables particle-resolved aerosol simulations using batched, scalable multiphase chemistry (Sect. 2.12). The following subsections provide implementation details for each component of the implementation.

Figure 1 Diagram of TChem-atm highlighting the customizability of the chemical mechanism and interoperability with numerous host model aerosol treatments. Further customization is supported for the selection of various Kokkos-enabled solvers, which allow performance portable computing across numerous platforms. ¹ Carbon Bond 05: Yarwood et al. (2005); ² SIMPOL: Pankow and Asher (2008); ³ CRACMM: Pye et al. (2023).

2.1 TChem-atm Provenance

TChem-atm follows the software architecture of TChem v2.0, sharing its Kokkos-based parallel model, CMake-based build system, and mechanism-independent kinetic-model abstraction. The analytic and numerical Jacobian capabilities in TChem-atm are adapted directly from TChem and rely on TINES/SUNDIALS for time integration and Newton/Krylov linear algebra. What is new in TChem-atm compared to TChem is the multiphase extension: CAMP-style gas-aerosol partitioning are added, new data structures (AerosolModelData and AerosolModelConstData) are introduced, and the RHS kernels now account for gas–aerosol mass transfer. These atmospheric-chemistry extensions are not present in TChem and required new mechanism parsing, new source-term kernels, and a generalized notion of “computational particles” to support sectional, modal, and particle-resolved host models.

TChem-atm does not depend on TChem as an external library. All numerical functionality required for atmospheric chemistry – Jacobian construction, batched RHS evaluation, and stiff ODE integration – is now provided by the standalone TINES library, which was separated from TChem expressly to support multiple projects. Only a few legacy utility routines from TChem remain in TChem-atm; the combustion-specific components of TChem are not needed. This separation allows TChem-atm to evolve independently while sharing a common numerical backbone with TChem through TINES.

TChem-atm incorporates the multiphase-chemistry abstractions originally introduced in CAMP, but it does not embed CAMP directly nor merge the two code bases. CAMP provides a mechanism-agnostic representation of gas–aerosol processes (e.g., partitioning, heterogeneous chemistry, condensed-phase reactions) and a unifying “computational particle” interface that allows chemistry to operate consistently across modal, sectional, and particle-resolved host models. TChem-atm adopts these abstractions and reimplements the corresponding source-term formulations using its own Kokkos-based backend.

The specific components inherited conceptually from CAMP include: (i) the representation of aerosol–gas coupling processes (mass transfer, heterogeneous reactions, etc.); (ii) the use of mechanism files to define multiphase chemistry at runtime; and (iii) the computational-particle abstraction that enables TChem-atm to work with different aerosol frameworks.

What is new in TChem-atm is the implementation of these concepts within a fully performance-portable infrastructure, including: (i) new data structures for aerosol-phase species and their device-resident constant data (AerosolModelData/AerosolModelConstData); (ii) GPU-accelerated source-term kernels for gas–aerosol coupling; (iii) integration with TINES and Kokkos for batched multiphase ODE solves; and (iv) a unified atmospheric-chemistry mechanism parser that extends TChem's gas-phase parser to support multiphase processes.

Thus, TChem-atm should be viewed not as a code-level merge of CAMP and TChem, but as a new atmospheric-chemistry engine that inherits CAMP's abstractions conceptually while implementing them natively within the TChem/TINES performance-portable numerical framework.

Table 1 Relationship between TChem, CAMP, and TChem-atm. Components are categorized as inherited as “Code” or “Concept”.

2.2 Chemical Mechanism Configuration

To support flexibility and ease of use, both CAMP and TChem-atm adopt human-readable configuration systems based on JSON or YAML for defining chemical mechanisms and model parameters (Dawson et al., 2022; Kim et al., 2023). This design enables users to specify complex reaction mechanisms, thermodynamic properties, and aerosol interaction parameters at runtime, without modifying the core codebase. Figure 2 shows an example gas-phase mechanism in YAML format for a simple system involving one reaction and three species ($A+B\to C$). This configuration is used to simulate the same chemical system independently across N computational cells, where each cell represents a distinct instance of the ODE system such as grid points in a host model or batches in a unit test.

The structured and editable format of these configuration files, along with the availability of validation tools, makes them accessible to a broad range of users – including those without specialized programming experience. In this work, we retain CAMP's approach for specifying multiphase chemistry processes such as gas–aerosol partitioning and condensed-phase reactions, and integrate it into the TChem-atm infrastructure. This unification allows researchers to flexibly configure chemically detailed multiphase systems in a portable, performance-optimized environment suited for both comprehensive atmospheric models and standalone simulations (e.g., chamber or flow-tube studies).

Figure 2 Example YAML input specifying a simple gas-phase mechanism consisting of a single reaction with three species ($A+B\to C$) used in simulations with N independent computational cells, each solving a separate instance of the ODE system.

2.3 Aerosol and Gas Kinetic Model Data

TChem-atm organizes the reaction mechanism information for gas and aerosol phases into two structured data objects as shown in Fig. 3. The gas phase data is stored in the Kinetic Model Constant Data (kmcd) object, while the aerosol phase data is stored in the Aerosol Model Constant Data (amcd) object.

Constructing these data objects involves two steps. First, the reaction mechanism is parsed and the relevant kinetic or thermodynamic information is stored in intermediate Kinetic Model and Aerosol Model data objects. These data are stored in dual Kokkos views, which allow synchronization of data between host and device memory. The data are further categorized and organized to optimize device-side kernel computation and facilitate postprocessing. For instance, for reactions following the Arrhenius form, TChem-atm stores both the reaction type and the associated parameters in separate dual views, as presented in Fig. 3.

Second, TChem-atm extracts only the data required by the computational kernels and synchronizes this subset to the device. The resulting constant data object – either kmcd or amcd – is then passed to the computation kernel. As a constant object, it cannot be modified during computation.

Figure 3 TChem-atm's internal representation of the chemical mechanism is the kinetic model object.

2.4 System state vector and governing equations

We define the system state vector η as the concatenation of the gas-phase state vector g and aerosol-phase chemical state vector μ such that $\mathit{\eta }=[\mathit{g},\mathit{\mu }]$. The gas-phase state is represented by mixing ratios of K chemical species $\mathit{g}=[g_{\mathrm{1}},g_{\mathrm{2}},\mathrm{\dots },g_K]$, while the aerosol phase is described by a set of N aerosol-phase vectors $\mathit{\mu }=[{\mathit{\mu }}_{\mathrm{1}},{\mathit{\mu }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\mu }}_N]$. Each aerosol-phase vector μ_i can be thought of as a “computational particle” i and contains the mass of A tracked species, ${\mathit{\mu }}_i=[{\mathit{\mu }}_i^{\mathrm{1}},{\mathit{\mu }}_i^{\mathrm{2}},\mathrm{\dots },{\mathit{\mu }}_i^A]$.

In this section, we assume a particle-resolved context for illustrative purposes, where each μ_i corresponds to a single particle. However, TChem-atm can also be applied to sectional and modal models, where the computational particles are then mapped to sections and representative particles in a mode. A detailed discussion of this mapping is provided in Sect. 2.7.

In expanded form, the system state vector η is

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{\eta }& =[\underset{\text{gas}}{\underbrace{g_{\mathrm{1}},g_{\mathrm{2}},\mathrm{\dots },g_K}},\underset{\text{particle 1}}{\underbrace{{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{1}},{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },{\mathit{\mu }}_{\mathrm{1}}^A}},\underset{\text{particle 2}}{\underbrace{{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{1}},{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}},\mathrm{\dots },{\mathit{\mu }}_{\mathrm{2}}^A}},\\ & \mathrm{\dots },\underset{\text{particle }N}{\underbrace{{\mathit{\mu }}_N^{\mathrm{1}},{\mathit{\mu }}_N^{\mathrm{2}},\mathrm{\dots },{\mathit{\mu }}_N^A}}].\end{array}\end{array}$$

In addition to the chemical state vector η, TChem-atm also requires the number concentration associated with each computational particle, denoted n_i for particle i. This number concentration is not part of the state vector because it is not modified within TChem-atm. Instead, it is passed in by the host model and may evolve independently due to physical processes such as coagulation, nucleation, or sedimentation. Nonetheless, n_i plays a critical role in the calculation of aerosol chemical source terms and is required for constructing the coupled gas–aerosol source terms described below.

The governing equations for the gas-phase and aerosol-phase components of the state vector are given by:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{g}_k}{\mathrm{d}t}}=F_{\mathrm{gg},k}+F_{\mathrm{ga},k},& {\displaystyle }k=\mathrm{1},\mathrm{\dots },K\text{(3)}& {\displaystyle \frac{\mathrm{d}{\mathit{\mu }}_i^a}{\mathrm{d}t}}=F_{\mathrm{ag},i}^a,& {\displaystyle }i=\mathrm{1},\mathrm{\dots },N,a=\mathrm{1},\mathrm{\dots },A\end{array}$$

The source term F_gg,k accounts for gas-phase chemical reactions that produce or consume species k, as specified by the gas-phase chemical mechanism and further explained in Sect. 2.5. The term F_ga,k represents gas–aerosol coupling processes, such as condensation of low-volatility species, evaporation of semi-volatile compounds, or uptake via heterogeneous reactions, and thus transfers mass between the gas phase and the aerosol phase. The corresponding aerosol-phase term, $F_{\mathrm{ag},i}^a$, describes the gain or loss of species a in computational particle i due to these interactions. These coupling terms ensure mass conservation across phases and encode the dynamic exchange of chemical species between gas and aerosol reservoirs and are explained in detail in Sect. 2.6. TChem-atm integrates the governing equations forward in time using the TINES or SUNDIALS CVODE ODE solver.

2.5 Gas Chemical Mechanism Source Term Construction

The gas-phase chemical source term, F_gg,k in Eq. (2), is computed via:

$$\begin{array}{}\text{(4)}& F_{\text{gg},k}=\sum _{j=\mathrm{1}}^{N_{\text{react}}}{\mathit{\nu }}_{kj}{q}_j,{\mathit{\nu }}_{kj}={\mathit{\nu }}_{kj}^{\prime \prime }-{\mathit{\nu }}_{kj}^{\prime },\end{array}$$

Here q_j is the reaction rate of reaction j, ${\mathit{\nu }}_{kj}^{\prime }$, ${\mathit{\nu }}_{kj}^{\prime \prime }$ are stoichiometric coefficients for species k in reaction j, and N_react is the total number of gas-phase reactions. The reaction rate is given by:

$$\begin{array}{}\text{(5)}& q_j=k_{\mathrm{f},j}\prod _{k=\mathrm{1}}^K{g}_k^{{\mathit{\nu }}_{kj}^{\prime }},\end{array}$$

where k_f,j is the reaction rate constant for reaction j. The functional forms for k_f,j depend on reaction type, as listed in Table 2.

Currently, TChem-atm can reproduce gas chemistry for two complex reaction mechanisms: the gas chemistry of E3SM v3, i.e., the UCI chemistry system (University of California Irvine), and the Carbon Bond 2005 chemical mechanism, which is well-formulated for urban to remote tropospheric conditions (Yarwood et al., 2005). Because TChem-atm is mechanism-agnostic, it can be easily adapted to support a wide range of applications. For example, a minimal chemical mechanism can be defined for instructional use in classroom settings, while detailed mechanisms with hundreds of species and reactions can be incorporated for high-fidelity research simulations.

Table 2 Gas reaction types available in TChem-atm; further details can be found in the TChem-atm online documentation.

2.6 Aerosol Mechanism Source Term Construction

TChem-atm computes the aerosol-gas interaction source terms, F_ga,k and $F_{\text{ag},i}^a$, as follows. In this section, we focus on the partitioning of semi-volatile gases. For simplicity, we assume that the first N_proc entries in the gas-phase species list correspond to partitioning species, and that each of these has a matching aerosol-phase counterpart with the same index a. That is, gas species g_a partitions to aerosol species ${\mathit{\mu }}_i^a$ for $a=\mathrm{1},\mathrm{\dots },N_{\mathrm{proc}}$.

$$\begin{array}{}\text{(6)}& {\displaystyle }{F}_{\text{ga},a}& {\displaystyle }=-\sum _{i=\mathrm{1}}^N{k}_i^a\left(g_a-g_a^{*}\right),& {\displaystyle }a& {\displaystyle }=\mathrm{1},\mathrm{\dots },N_{\mathrm{proc}},\text{(7)}& {\displaystyle }{F}_{\text{ga},k}& {\displaystyle }=\mathrm{0},& {\displaystyle }k& {\displaystyle }=N_{\mathrm{proc}}+\mathrm{1},\mathrm{\dots },K,\text{(8)}& {\displaystyle }{F}_{\text{ag},i}^a& {\displaystyle }={\displaystyle \frac{{\mathit{\rho }}^a}{n_i}}{k}_i^a\left(g_a-g_a^{*}\right),\end{array}$$

where $k_i^a$ (s⁻¹) is the mass transfer coefficient for species a and particle i, $g_a^{*}$ (mole mole⁻¹) is the equilibrium mixing ratio of species a, ρ^a (kg m⁻³) is the mass density of gas species a needed for the conversion of mixing ratio to mass concentration, and n_i is the number concentration of computational particle i needed to convert mass concentration to mass. While F_ga,a is the mixing ratio flux of species a to the gas phase from all particles, $F_{\text{ag},i}^a$ is the mass flux of species a to computational particle i.

The mass transfer coefficient $k_i^a$ and the mass density ρ^a are given by,

$$\begin{array}{}\text{(9)}& {\displaystyle }{k}_i^a& {\displaystyle }=\mathrm{4}\mathit{\pi }{r}_i{D}_{\mathrm{g}}^a{n}_{i}f(K{n}_i^a,{\mathit{\alpha }}_a),\text{(10)}& {\displaystyle }{\mathit{\rho }}^a& {\displaystyle }={\displaystyle \frac{p{M}^a}{RT}},\end{array}$$

where r_i is the effective radius of the particle, $D_{\mathrm{g}}^a$ is the gas-phase diffusion coefficient of species a, and $f(K{n}_i^a,{\mathit{\alpha }}_a)$ is the Fuchs-Sutugin correction factor, p is atmospheric pressure, M^a is the molar weight of gas a (kg mole⁻¹), R is the universal gas constant (J mole⁻¹ K⁻¹), and T is temperature (K).

The effective radius r_i can be derived from the particle masses ${\mathit{\mu }}_i^a$ and aerosol species densities ${\mathit{\rho }}_{\mathrm{p}}^a$ by assuming sphericity:

$$\begin{array}{}\text{(11)}& r_i=\left({\displaystyle \frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}}\sum _{a=\mathrm{1}}^A{\displaystyle \frac{{\mathit{\mu }}_i^a}{{\mathit{\rho }}_{\mathrm{p}}^a}}\right).\end{array}$$

The Fuchs-Sutugin correction factor is computed as:

$$\begin{array}{}\text{(12)}& f_{\mathrm{fs}}(K{n}_i^a,{\mathit{\alpha }}_a)={\displaystyle \frac{\mathrm{0.75}{\mathit{\alpha }}_a(\mathrm{1}+K{n}_i^a)}{(K{n}_i^a{)}^{\mathrm{2}}+(\mathrm{1}+\mathrm{0.283}{\mathit{\alpha }}_a)K{n}_i^a+\mathrm{0.75}{\mathit{\alpha }}_a}},\end{array}$$

where $K{n}_i^a={\mathit{\lambda }}_a/r_i$ is the Knudsen number, λ_a is the mean free path of gas a, and α_a is the mass accommodation coefficient of gas a, here assumed to be 0.1 for all gases.

To calculate the equilibrium mixing ratios g^* we use the SIMPOL.1 method (Pankow and Asher, 2008), which parameterizes $g_a^{*}$ as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{g}_a^{*}& {\displaystyle }={\displaystyle \frac{p_a^{*}}{p}},\text{(14)}& {\displaystyle }{\log }_{\mathrm{10}}{\displaystyle \frac{p_a^{*}}{p_{\mathrm{0}}}}& {\displaystyle }={\displaystyle \frac{B_{\mathrm{1}}^a}{T}}+B_{\mathrm{2}}^a+B_{\mathrm{3}}^{a}T+B_{\mathrm{4}}^a\ln \left(T\right),\end{array}$$

where $p_a^{*}$ is the equilibrium vapor pressure for gas a (Pa), p₀=101 325.0 Pa is standard atmospheric pressure and $B_{\mathrm{1}}^a$, $B_{\mathrm{2}}^a$, $B_{\mathrm{3}}^a$, $B_{\mathrm{4}}^a$ are fitting parameters provided by SIMPOL.1.

This section focuses on semi-volatile partitioning because this is the only heterogeneous process type currently implemented in TChem-atm. The framework, however, is extensible: heterogeneous reactions, reactive uptake parameterizations, or condensed-phase chemical processes can be incorporated by defining the corresponding rate expressions in the mechanism file and adding the associated source-term kernels. These capabilities follow the CAMP abstraction and will be implemented in future extensions of TChem-atm. At present, semi-volatile gas–aerosol partitioning represents the primary heterogeneous pathway supported.

2.7 Abstracting Aerosol Representations Across Models

One of CAMP's key innovations is to perform chemistry on individual computational particles, providing a unified library that works across different aerosol representations. As previously explained, a computational particle may represent an actual particle (in particle-resolved models), a bin section (in sectional models), or a representative particle within a mode (in modal models). This approach replaces the traditional coupling of chemistry to specific aerosol schemes – which is often hard-coded and inflexible – with a more general method that simplifies the addition of new species or processes. While this idea was originally described as an “abstraction of aerosol representation” in Dawson et al. (2022), the current framing in terms of computational particles offers a more intuitive and practical perspective, and steps toward the generalized representations envisioned in Hodzic et al. (2023).

Figure 4 Representation of the mapping between the abstract representation of the aerosol state vector μ in TChem-atm and the host model representation (blue filled regions). The host model is responsible for tracking the number concentration of each computational particle, contained in the n vector. Particle diameters are diagnosed from μ and n.

Following this approach, TChem-atm defines a general interface for aerosol–gas chemical interactions. Rather than prescribing how aerosols must be represented, TChem-atm relies on the host model to provide two essential pieces of information for each computational particle: (1) the chemical composition (i.e., the vector of aerosol species masses), and (2) the associated number concentration. As noted in Sect. 2.4, the number concentration is not part of the TChem-atm state vector, since it is not modified by the chemistry, but it plays a critical role in constructing source terms and must be supplied by the host model. Similarly, the particle diameter is required for processes such as surface-area-based uptake and size-dependent reaction rates. It is diagnosed from the particle masses and the number concentrations (assuming spherical particles). In particle-resolved models, this diameter is simply the diameter of the individual computational particle. In sectional models, it is typically taken as the midpoint of the bin.

For modal models, however, the situation is more complex. Modal schemes typically do not pass explicit diameters to the chemistry module; instead, they compute the effects of condensation by integrating over assumed size distributions and updating moments analytically – an approach rooted in the work of Whitby and McMurry (1997). While this method is widely used, it complicates the integration with approaches like TChem-atm that expect a well-defined particle diameter for each computational particle. This limitation is not a fundamental barrier, but addressing it would require additional work. Specifically, modal host models could be adapted to pass representative diameters (e.g., mean or characteristic values), which would improve compatibility with TChem-atm and enhance transparency and efficiency in multiphase chemistry calculations.

TChem-atm requires that each computational particle contains the full set of aerosol-phase species defined by the chemical mechanism. While individual species masses may be zero in any given computational particle, all species must be present in the data structure. This design ensures a consistent state vector structure across computational particles and simplifies source term construction. As a result, TChem-atm does not currently support mode-specific species sets, such as those used in some modal models where particular modes are restricted to subsets of species (e.g., black carbon only, or sea salt with organic carbon only). Any such mapping or aggregation must be handled externally by the host model prior to passing data to TChem-atm.

The abstractions within TChem-atm also facilitate integration with emerging or hybrid aerosol representations, such as adaptive sectional models, quadrature-based methods, or machine-learning-driven representations, as long as the host model can supply the required species masses, number concentrations, and diameters for each computational particle. Figure 4 illustrates how computational particles in TChem-atm correspond to different host model representations. For sectional and particle-resolved methods, there is a one-to-one mapping between computational particles tracked by TChem-atm and the corresponding host representation. For the modal approach, we display a scenario in which each mode is represented by a set of three characteristic computational particles with diameters corresponding to the geometric mean diameter and one standard deviation above and below.

2.8 Performance Portability with Kokkos

To enable performance portability across diverse computing platforms, PartMC-TChem-atm leverages the Kokkos library (Edwards et al., 2014; Trott et al., 2021, 2022). Kokkos is a programming model designed to abstract parallel execution and memory management, allowing a single implementation of numerical algorithms to run efficiently on a variety of architectures, including multi-core CPUs (OpenMP), NVIDIA GPUs (CUDA), AMD GPUs (HIP), and Intel GPUs (SYCL). Kokkos achieves this through compile-time polymorphism and execution policies that separate algorithm logic from hardware-specific optimizations. This abstraction is critical for modern atmospheric modeling, where simulations must scale across heterogeneous high-performance computing (HPC) systems.

Kokkos has been adopted in several modeling efforts to modernize physics and chemistry components for exascale computing. It underpins performance-portable implementations in the Energy Exascale Earth System Model (E3SM) atmosphere component (Taylor et al., 2023) and serves as the parallel backend for TChem (Kim et al., 2023). By integrating TChem-atm with Kokkos, PartMC-TChem-atm inherits the ability to execute detailed multiphase chemistry on both CPU and GPU architectures without modifying core algorithms.

2.9 Software Architecture and Jacobian Computation

TChem-atm adopts the same software design as the TChem library (Kim et al., 2023) and uses batched computations and algorithm implementations based on two primary Kokkos abstractions, ValueType and DeviceType. ValueType represents built-in scalar data types such as int and float, but it can also use automatic differentiation (AD) types from, for example, the SACADO package to support automatic computation of Jacobians and to facilitate differentiability for machine learning. DeviceType defines a memory space (where data lives, e.g., on CPU or GPU) and an execution space (where kernels run, e.g., on CPU or GPU). Each Kokkos backend supports a specific DeviceType so that applications written in Kokkos, using these abstractions, can run on diverse architectures. The workflow of TChem-atm is shown in Algorithm 1. It creates two data objects for gas and aerosol chemistry and uses hierarchical parallelism to batch-evaluate source terms or initial value problems. Nested team-thread parallelism enables performance across varying system configurations.

TChem-atm automatically computes Jacobian matrices using either finite difference methods or automatic differentiation via SACADO (Phipps et al., 2020). When ValueType is double, a numerical Jacobian is computed; TINES supports forward differencing, central differencing, and Richardson extrapolation, with forward differencing used by default. For analytical Jacobians, ValueType is an automatic differentiation type, and SACADO computes derivatives using the chain rule. This unified design simplifies the interface for downstream solvers.

As previously discussed, TChem-atm leverages Kokkos' team parallelism, utilizing batched computation. Specifically, multiple evaluations of nested routines to invoke a team function. The number of batches or samples is directly correlated with the number of grid points, as TChem-atm is invoked by an atmospheric model, or with model evaluations, where TChem-atm is utilized for uncertainty quantification and machine learning studies. Within the batched loop, a team routine is executed, with the primary team routine in TChem-atm being the right-hand-side routine, as shown in Algorithm 1.

For the gas component of the source term, nested parallel loops are employed to compute the kinetic constants for each reaction type (see Table 2), the rate of progress for each reaction, and the net production rate. For the aerosol component, nested parallel loops are similarly utilized for the evaluation of each aerosol process, where the number of particles or sections determines the size of the loop; currently, the SIMPOL process (Pankow and Asher, 2008) is supported within TChem-atm. Due to this batched design, we do not anticipate superior computational performance from TChem-atm when only a single evaluation is conducted. However, enhanced performance is expected as the number of evaluations increases.

Algorithm 1 Batched algorithm for TChem-atm computation of chemical source terms.

2.10 ODE Solvers and Integration Workflow

TChem-atm is designed to decouple the specification of chemical processes from the numerical integration of the resulting ODE systems. In atmospheric models, gas- and aerosol-phase chemistry are typically solved independently at each grid point, often requiring the integration of stiff ODE systems over short time scales.

To support this need, TChem-atm provides interfaces to two widely used implicit solvers: TINES (Kim and Diaz-Ibarra, 2021) and SUNDIALS CVODE (Gardner et al., 2022; Hindmarsh et al., 2005; Balos et al., 2021; Hindmarsh et al., 2025), both of which implement backward differentiation formula (BDF) methods suitable for stiff systems. TINES employs a second-order Trapezoidal BDF (TrBDF2) scheme (Bank et al., 1985) with a Newton-based nonlinear solver and UTV decomposition for solving linear systems. CVODE supports both dense LU-based solvers, which require a Jacobian matrix and can use either analytical Jacobians (computed via automatic differentiation with SACADO (Phipps et al., 2020)) or numerical Jacobians (via finite difference schemes) (Salane, 1986), and GMRES-based iterative solvers, which do not require explicit Jacobians and instead approximate the action of the Jacobian through matrix-free methods.

To enable scalable integration across many independent ODE systems, TChem-atm adopts a batched interface, where source term evaluations and Jacobian computations are executed in parallel using Kokkos' team-level parallelism. This structure is compatible with large-scale simulations that solve thousands to millions of ODE systems simultaneously, such as grid-point chemical updates in 3D models.

Users can configure solver settings – including absolute and relative tolerances, maximum time steps, and method-specific options – via human-readable input files. This flexibility, along with modular interfaces to both solvers, allows users to tailor solver behavior to problem stiffness and computational constraints without modifying the underlying chemical mechanism or solver logic.

The interchangeable solver design also facilitates future extensibility. Additional solvers can be integrated with minimal code changes, supporting experimentation with different numerical strategies or coupling to new host models. Together, these features make TChem-atm a flexible and performance-portable platform for integrating chemically detailed, stiff ODE systems in a wide range of atmospheric chemistry applications. Performance comparisons across solver configurations – including dense and Jacobian-free methods on CPU and GPU architectures – are discussed in Sect. 3.2.

2.11 Testing Strategy

Development of TChem-atm is maintained via a GitHub repository. To ensure changes to the code base do not break model execution, numerous tests are run via GitHub Actions (GitHub, 2025) upon both the opening and merging of pull requests. Tests include small-scale unit testing of reaction and production rates, as well as larger-scale integration tests in which the entire chemical mechanism including E3SM's UCI chemistry system and the Carbon Bond 05 mechanism are evaluated. Additionally, construction of the source term is tested for the SIMPOL gas-aerosol system (Pankow and Asher, 2008) which describes partitioning of organic species.

2.12 Enabling Particle-Resolved GPU Chemistry

We integrated TChem-atm into the PartMC particle-resolved model to enable GPU-accelerated, chemically detailed simulations of multiphase aerosol chemistry. PartMC is a stochastic, particle-resolved aerosol box model that resolves the composition of many individual aerosol particles within a well-mixed volume of air. Riemer et al. (2009), DeVille et al. (2011), Curtis et al. (2016), and DeVille et al. (2019) describe in detail the numerical methods used in PartMC. PartMC tracks a large number of discrete computational particles to represent the aerosol population. Each computational particle is a vector of composition, which tracks the mass of each species within the particle and evolves the composition in time. By coupling PartMC with TChem-atm, each particle's chemistry is advanced using TChem-atm's batched source term evaluations and solver interface. Kokkos ensures that these evaluations are offloaded to GPUs with minimal overhead.

This integration allows PartMC to take advantage of the performance portability and multiphase chemistry capabilities developed in this work. The interface supports flexible mapping between PartMC's particle state variables and TChem-atm's kinetic model structure, preserving modularity and extensibility. In test cases, we validated the consistency of results across CPU and GPU architectures.

3 Results

This section presents two sets of results. All simulations underlying the results presented here are fully reproducible using the archived code and dataset described in Diaz-Ibarra et al. (2025a, b, c). Because the goal of this work is to establish and validate the TChem-atm framework, the tests presented here focus on controlled box-model scenarios that isolate chemical behavior and computational scaling. First, we evaluate the chemical accuracy of the new TChem-atm integration within PartMC by comparing its outputs against the established PartMC-CAMP model in a box model simulation (Sect. 3.1). This includes gas-phase chemistry and gas–aerosol partitioning for a previously publised scenario (Dawson et al., 2022). Second, we assess the computational performance of the TChem-based implementation across multiple hardware backends (CPUs and GPUs) and solver configurations. We report per-particle wall-clock timings for both full ODE solves and individual kernel components such as RHS and Jacobian evaluations. Together, these results demonstrate the correctness and scalability of PartMC–TChem-atm for atmospheric chemistry modeling.

For this study, PartMC-TChem-atm was compiled and tested using Kokkos 4.5.01 with support for CPUs and GPUs. Simulations were performed on the following computing architectures: an NVIDIA H100 GPU on DeltaAI at the National Center for Supercomputing Applications (NCSA), University of Illinois; an AMD MI250x GPU on Frontier at the Oak Ridge Leadership Computing Facility (OLCF); and an AMD EPYC 7763 CPU on Perlmutter at the National Energy Research Scientific Computing Center (NERSC). The build process used the CMake-based infrastructure provided by TChem-atm, with Kokkos execution and memory spaces selected at compile time to target the appropriate backend. GPU simulations used the CUDA and HIP backends, while CPU runs used the OpenMP backend. All simulations used identical chemical mechanisms, numerical tolerances, and time step settings to isolate hardware performance and consistency of results.

3.1 PartMC-TChem-atm box model setup and comparison with PartMC-CAMP

To evaluate the performance of TChem-atm, we conducted a box model simulation using PartMC. The model was configured to represent a well-mixed air parcel with constant temperature and pressure. This simulation focuses exclusively on gas-phase chemistry and gas–aerosol partitioning, with coagulation and removal processes disabled. Gas-phase reactions were governed by the CB05 mechanism (Yarwood et al., 2005), and gas–aerosol partitioning followed the SIMPOL scheme (Pankow and Asher, 2008). The model was run for 24 h with a time step of 60 s. Initial gas-phase concentrations and emissions were set according to Table 10 of Dawson et al. (2022), while the aerosol phase was initialized with three modes, as described in the same reference. For the simulations presented here, we initialized PartMC with 10³ computational particles. The complete set of code, model configuration files, initial conditions, and post-processing scripts used for this experiment are archived and publicly available (Diaz-Ibarra et al., 2025a, b, c).

Simulation results are shown in Fig. 5. Panel (a) displays the 24 h time series of ozone mixing ratios as predicted by PartMC-CAMP, PartMC-TChem-atm (CPU), and PartMC-TChem-atm (GPU). Panels (b)–(d) shows the evolution of gas-phase isoprene (ISOP) and its semi-volatile products, ISOP-P1 and ISOP-P2. ISOP-P1 is formed through the reaction of ISOP with OH, while ISOP-P2 results from its reaction with ozone. Figure 5e–f presents the aerosol-phase mass concentrations of ISOP-P1_aero and ISOP-P2_aero, which arise from the partitioning of their respective gas-phase precursors.

Figure 5 Comparison between PartMC-TChem-atm (CPU), PartMC-TChem-atm (GPU) and PartMC-CAMP for (a) ozone mixing ratio, (b) ISOP, (c) ISOP-P1 and (d) ISOP-P2 mixing ratios, (e) ISOP-P1_aero and (f) ISOP-P2_aero mass concentration for the 24 h simulation period.

Per-timestep relative tolerances provided by TChem-atm to the CVODE solver control the final accumulated error. In order to achieve 0.01 % global accuracy for this comparison, a relative tolerance of 10⁻⁶ was chosen. To summarize how PartMC-TChem-atm compares to PartMC-CAMP, the ℓ²-norm of the relative difference was $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}}$ for gas phase species and $\mathrm{8.2}\times {\mathrm{10}}^{-\mathrm{6}}$ for aerosol species.

To more clearly quantify the agreement between the CPU, GPU, and legacy PartMC–CAMP implementations, Fig. 6 shows the relative differences between the CPU, GPU, and PartMC–CAMP results. The discrepancies remain small throughout the simulation, with a modest temporal drift for some species that reflects normal numerical variability rather than any architecture-specific effect. All deviations remain well within solver tolerances.

Overall, PartMC-TChem-atm reproduces the results of PartMC-CAMP with high fidelity. It performs consistently across both CPU and GPU architectures, capturing expected chemical behavior while enabling accelerated simulations on heterogeneous computing platforms.

Figure 6 Relative differences between CPU, GPU, and PartMC–CAMP mixing ratios for the species shown in Fig. 5.

3.2 Benchmarking GPU Speedup and Solver Scaling

We evaluate the computational performance of TChem-atm on several hardware backends including multicore CPUs and both NVIDIA and AMD GPU architectures. To achieve this, the box model example from Sect. 3.1 was run by varying the number of computational particles, which increases the size of the source term and the Jacobian matrix. We tested several ODE solver configurations presented in Table 3, including TINES TrBDF2 and two variants of SUNDIALS CVODE that differ by the sparsity of the linear solver (SUNDIALS CVODE-GMRES and SUNDIALS CVODE (dense)). The same numerical parameters are applied for the three configurations.

Figure 7 Per-particle solver wall clock time through t=10 s organized by solver. Timings are colored by computing architecture. For each solver and architecture, team size and vector size were set via the Kokkos::AUTO configuration. TINES TrBDF2 timings for the AMD MI250X at particle counts larger than N_p∼100 are not included due to the required memory exceeding allocatable GPU DRAM capacity. Shaded regions indicate approximate computing regimes for each aerosol treatment.

For all numerical experiments, the model was run to a simulated time of 10 s and the wall clock time per time step was measured. For each solver, the number of computational particles was varied from N_p=1 to $N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ in twenty, logarithmically spaced increments. Figure 7 shows the wall clock time per particle for each solver configuration across a range of particle counts. For small particle numbers (i.e., N_p≲50), solvers using dense linear algebra – namely SUNDIALS CVODE (dense) and TINES TrBDF2 – achieve the best performance due to their low overhead and efficient handling of small Jacobians. However, as N_p increases, the cost of forming and solving dense Jacobian systems grows rapidly, making these approaches computationally prohibitive beyond N_p∼500.

In contrast, SUNDIALS CVODE-GMRES employs a Jacobian-free Krylov method, which avoids explicit construction of the Jacobian. This enables continued scaling to larger particle numbers, with wall clock time per particle decreasing even up to N_p=10⁶ on GPU architectures. These results underscore the importance of matching solver structure to problem size and hardware – dense solvers excel at small scale, while Jacobian-free methods offer superior scalability.

In the GMRES case, we observe substantial speedups of GPU relative to CPU at large computational particle numbers, underscoring the advantages of Jacobian-free Krylov methods on GPUs. These results are based on a single-cell configuration, however, and therefore do not yet fully exploit the parallelism available on GPU architectures. In larger multi-cell simulations, we anticipate even greater gains for GMRES as well as improvements for dense solvers, which are not fully reflected under the present single-cell conditions.

Table 3 Numerical parameters for evaluated ODE solvers.

These results have practical implications for model selection and hardware deployment. Particle-resolved simulations such as PartMC typically require on the order of 10 000 computational particles per grid box to adequately resolve mixing state and aerosol composition. In this regime, dense linear solvers become inefficient, and the GPU implementation of SUNDIALS CVODE-GMRES offers significant performance advantages due to its Jacobian-free formulation and better scalability with particle number.

Conversely, if the application involves simplified aerosol representations – such as modal or sectional models with a small number of size/composition bins – the computational burden is much lower, and CPU-based solvers with dense linear algebra may offer better performance. Thus, the optimal hardware–solver combination depends strongly on the complexity and resolution of the aerosol representation. To further illustrate the relationship between aerosol treatment and the choice of solver, approximate computing regimes for each aerosol treatment are shown as shaded regions in Fig. 7.

Timings shown in Fig. 7 were conducted using the Kokkos:AUTO configuration, which automatically sets the level of parallelism at both the team and vector level (for those familiar with NVIDIA systems, the Kokkos team size and vector size are cast as the x and y dimension of each CUDA block). Both parameters can be set at runtime to determine a problem-specific optimal configuration for team and vector size. Figure 8 shows heatmaps of wall clock time per particle for the SUNDIALS CVODE-GMRES solver where both the team size and vector size are varied. Optimal team and vector size configurations differ by GPU. For the NVIDIA H100, a team size of 16 and a vector length of 32 resulted in the lowest wall clock time per particle, whereas a team size of 64 and vector length of 16 were best for the AMD MI250X. Below each heatmap in Fig. 8 are the minimum, maximum, and mean wall clock time across each team and vector configuration alongside the wall clock time achieved by the Kokkos:AUTO configuration. We find that the optimal configuration for the NVIDIA H100 is within a factor of 2 faster than the Kokkos:AUTO configuration, while both are equal for the AMD MI250X.

Figure 8 Heatmaps of wall clock time per particle as a function of Kokkos team size and vector length for the SUNDIALS CVODE-GMRES solver on the NVIDIA H100 (left) and AMD MI250X (right). For each solver run, the number of particles was fixed at 1000. The optimal configuration (lowest wall clock time) is indicated by a white star in each plot.

Ultimately, solver performance is governed by the size and structure of the Jacobian matrix. A large number of particles or reactions leads to a large Jacobian, while a large number of chemical species often results in a sparse one. The balance between these factors – reaction count, species count, and particle/bin resolution – ultimately determines whether dense or iterative solvers, and CPU or GPU hardware, offer the most efficient path forward.

In addition to full ODE solver runs, a separate set of timings were conducted to evaluate the performance of the RHS evaluation and Jacobian construction. Figure 9 shows how wall clock time for the RHS and Jacobian vary as the number of particles increases.

For each computing backend, per-particle performance for RHS evaluations improves as the number of particles increases. When $N_{\mathrm{p}}\lesssim \mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$, we find the best performance for the NVIDIA H100. Above $N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$, the AMD MI250X achieves the best performance through $N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$. At high particle counts ($N_{\mathrm{p}}\gtrsim \mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}$), performance plateaus for each computing backend. At $N_{\mathrm{p}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$, GPU backends achieve speedups of 2.3 and 5 relative to the CPU for the NVIDIA H100 and AMD MI250X, respectively.

Additionally, Jacobian evaluation timings exhibit similar performance for both the CPU and GPUs: wall clock time per particle decreases when the number of particles is low (i.e., N_p≲100), however, performance degrades as N_p increases further due to corresponding growth in the Jacobian.

To interpret these performance trends, it is helpful to clarify the role of direct versus iterative solvers. We include dense-direct linear solvers in our performance evaluation not because they are expected to be optimal for large atmospheric chemistry systems – indeed, the Jacobians arising from CB05 and similar mechanisms are sparse – but because they provide a well-defined performance baseline for small to moderate system sizes. Dense solvers incur higher memory costs and scale poorly as the number of particles grows, so their performance degradation is expected. However, these solvers remain competitive in regimes with few computational particles or simplified aerosol representations (e.g., small sectional or modal systems), where the Jacobian dimension remains modest and overhead is low. Including them therefore helps delineate the transition point at which preconditioned iterative methods such as GMRES become the clearly superior choice. Our intention is not to advocate dense solvers for large, stiff atmospheric chemistry problems, but rather to provide a complete performance map across solver categories and hardware backends. For realistically sized atmospheric mechanisms, preconditioned iterative Krylov methods are preferable because they exploit sparsity and dramatically reduce memory demands. The results in Sect. 3.2 are therefore consistent with expectations and are included to contextualize solver performance across operational regimes.

Figure 9 Per-particle wall clock time for evaluation of the gas-aerosol right-hand side (left) and Jacobian (right). Timings are colored by computing architecture. For each timing run, team size and vector size were set via the Kokkos::AUTO configuration. Particle number is limited for Jacobian evaluations due to exceeding the allocatable memory of each computing architecture.

We also conducted lightweight profiling during the timing experiments to separate kernel execution from host–device memory transfer. For the idealized box-model runs presented here, more than 99.9 % of the wall-clock time occurs inside the Kokkos kernels, with memory transfers accounting for only the remaining fraction. This behavior is expected in the present setup, where particles share identical initial conditions, the RHS evaluation is uniform across particles, and memory movement is minimized. In more realistic host-model applications, such as full PartMC–TChem-atm simulations, additional data exchange between the host model and TChem-atm, more heterogeneous particle states, and increased stiffness may increase the relative cost of memory transfers and solver iterations. Thus, while our profiling confirms that kernel execution dominates in these proof-of-concept tests, the performance breakdown will necessarily depend on the host-model coupling pattern and the chemical and microphysical variability of the simulated system.

3.3 Potential Applications in Chemistry–Transport and Climate Models

While the tests presented in this paper focus on idealized box-model configurations to isolate numerical behavior and performance, TChem-atm is designed for use in full chemistry–transport and Earth system models. CAMP, from which TChem-atm adopts its multiphase abstractions, has already been integrated into the MONARCH model system (Dawson et al., 2022), demonstrating the suitability of the CAMP structure for real-world chemical forecasting and case-study applications. To date, TChem-atm has not yet been deployed in operational chemistry–transport models; however, the TChem-atm interface is compatible with the operator-splitting approach used in such models, where chemistry is solved independently in each grid cell as a box model. Integration therefore requires only the development of a host-model-specific interface, not changes to TChem-atm itself.

A near-term target for deployment is the Energy Exascale Earth System Model (E3SM), where the Kokkos-based design of TChem-atm is well aligned with ongoing efforts to modernize the atmospheric chemistry infrastructure for heterogeneous computing environments. Exploring this integration constitutes part of our future work and leverages the fact that TChem-atm already implements the data structures and parallelism patterns required for grid-level chemistry calls in a chemistry–transport or climate model.

4 Conclusions

We have presented TChem-atm, a performance-portable approach for simulating chemically detailed, multiphase atmospheric chemistry across diverse computing architectures. By integrating the chemical mechanism infrastructure of CAMP with the Kokkos-based backend of TChem, TChem-atm enables consistent deployment across CPUs and GPUs while maintaining flexibility in mechanism configuration and solver selection.

A key strength of TChem-atm is its support for multiple aerosol representations – including modal, sectional, and particle-resolved approaches – making it suitable for a wide range of host models and applications. Integration with the particle-resolved model PartMC demonstrates that TChem-atm accurately reproduces gas-phase chemistry and gas–aerosol partitioning results, verifying the method's correctness and interoperability.

Performance benchmarks reveal that dense linear solvers are well-suited for small particle populations, whereas Jacobian-free solvers such as CVODE-GMRES offer scalable performance benefits on GPUs as the number of computational particles increases. These gains are particularly significant for particle-resolved models like PartMC, which typically require around 10 000 particles per grid cell to capture aerosol mixing state. Across solvers, our modular approach supports efficient, scalable integration of detailed multiphase chemistry across diverse computing platforms.

While the current version of TChem-atm is not yet fully performance-optimized, it provides a solid and extensible foundation for future development. Ongoing efforts will focus on enhanced solver strategies, improved memory layout, and broader chemical mechanism support. The code and all associated test cases are publicly available under an open-source license, enabling verification, reuse, and extension by the broader atmospheric modeling community (Diaz-Ibarra et al., 2025a, b, c).