Isentropic mixing vs. convection in CLaMS-3.0/MESSy: evaluation using satellite climatologies and in situ carbon monoxide observations

Lagrangian modeling of transport, as implemented in the Chemical Lagrangian Model of the Stratosphere (CLaMS), connects the advective (reversible) component of transport along 3D trajectories with mixing, the irreversible component. Here, we investigate the interplay between strongly localized convective uplifts and large-scale flow dynamics in the upper troposphere and lower stratosphere (UTLS). We revisit the Lagrangian formulation of convection in CLaMS-3.0/MESSy, driven by ECMWF's ERA5 reanalysis, and further develop the model. These developments include refining spatial resolution in the Planetary Boundary Layer (PBL) and decoupling the frequency of the adaptive grid procedure – which captures isentropic mixing and redefines Lagrangian air parcels – from the parameterization of convection.

To improve the model's UTLS transport representation, particularly from the PBL over days to weeks, we derive zonally and seasonally resolved climatologies of CO partial columns (XCO, spanning 147–68 hPa) and compare them with Microwave Limb Sounder (MLS) and Atmospheric Chemistry Experiment Fourier Transform Spectrometer (ACE-FTS) observations, as well as in situ data. Incorporating a parameterization for unresolved convection significantly improves CO anomaly representation in the UTLS, particularly in capturing seasonal and spatial patterns. While the simulated absolute XCO values align better with ACE-FTS, the model reproduces MLS anomalies more accurately, suggesting MLS better represents CO variability. In situ observations in the boreal polar region generally support lower ACE-FTS CO values, while MLS better represents CO enhancements in air affected by the Asian summer monsoon above 10 km.

1 Introduction

The Chemical Lagrangian Model of the Stratosphere (CLaMS) is a Chemistry Transport Model that employs an irregular Lagrangian grid, trajectory-based 3D advection driven by diabatic, cross-isentropic vertical velocities, and an adaptive grid procedure to simulate isentropic mixing. In this study, CLaMS is driven by ERA5 reanalysis data from ECMWF (Hersbach et al., 2020). Initially formulated as a purely stratospheric, isentropic model on selected potential temperature (θ) levels, early CLaMS applications focused on ozone depletion in the polar vortex (McKenna et al., 2002a, b).

Despite the clear advantages of the Lagrangian approach, which inherently separates mixing from the advective part of transport, very few global Lagrangian chemical transport models with explicit mixing exist (Reithmeier and Sausen, 2000; Stenke et al., 2008; Wohltmann and Rex, 2009; Pugh et al., 2012). One of the key challenges in atmospheric transport modeling is the accurate representation of mixing and convection, particularly in the troposphere and at the interface with the stratosphere. Recent studies have shown that the persistent moist bias in the lowermost stratosphere observed in current climate models (Stenke et al., 2008, 2009), often exceeding 200 %, is strongly influenced by the choice of transport schemes. Eulerian transport models tend to introduce excessive numerical diffusion, whereas less diffusive Lagrangian schemes, such as the one implemented in CLaMS, significantly reduce this bias, leading to improved climate simulations (Charlesworth et al., 2023; Ploeger et al., 2024). The latest version, CLaMS-2.0/MESSy, extends the model to include the troposphere and utilizes the Modular Earth Submodel System (MESSy) as its software framework (Jöckel et al., 2010). This version also runs on the JUWELS supercomputer (Jülich Supercomputing Centre, 2019; Konopka et al., 2019, 2022).

The aim of this paper is to revisit and improve the model formulation of unresolved convection, particularly its spatial and temporal resolution, and its interaction with the CLaMS mixing parameterization, which is based on the adaptive grid procedure (Konopka et al., 2004; Pommrich et al., 2014). In the first part of the paper, we review all CLaMS components related to transport processes. In the second part, we evaluate this new model version, CLaMS-3.0/MESSy, using a particularly well-suited diagnostic to assess transport in the Upper Troposphere and Lower Stratosphere (UTLS). Specifically, we compare partial columns of CO (spanning 147–68 hPa) observed by the Microwave Limb Sounder (MLS) and the Atmospheric Chemistry Experiment Fourier Transform Spectrometer (ACE-FTS), as well as in situ observations, with corresponding model data (Wright et al., 2025). Our analysis focuses on zonally resolved seasonal CO patterns, averaged over multiple years, which are largely influenced by seasonal variability in convection. Since chemistry plays a minor role in this study, we employ a simplified CLaMS configuration that includes CO, ozone, and water vapor, with boundary conditions as described by Pommrich et al. (2014).

2 CLaMS adaptive grid with convection

Unlike the more common Eulerian frame of reference, where atmospheric motion is observed from a fixed spatial position, the Lagrangian approach follows individual fluid parcels as they move with the flow. In a numerical model, this distinction means that while the Eulerian grid remains fixed in space, the Lagrangian grid moves along with the flow, materially anchored to it. Thus, the Lagrangian grid, or Lagrangian air parcels (APs), can be understood as massless pivotal points that are intrinsically connected with the flow. Their motion, driven by the flow velocities – either forward or backward in time – defines the so-called Lagrangian trajectories. In both frames of reference, the physical properties at a given grid point represent the mean properties of the atmosphere within the volume of the grid box surrounding that point. While calculating grid volume in Eulerian models is straightforward, Lagrangian grids require the concept of nearest neighbors and Voronoi volumes (see Appendix A).

Lagrangian transport is inherently less diffusive than Eulerian transport, the latter requiring frequent interpolations on the background grid due to the Courant–Friedrichs–Lewy condition (Courant et al., 1928), which often leads to substantial numerical diffusion, exceeding observed atmospheric diffusion. However, purely Lagrangian trajectories without any mixing represent an unrealistic extreme that does not align with observations (see e.g. Pan et al., 2006). CLaMS addresses this limitation through an adaptive grid procedure that combines Lagrangian trajectories (i.e. the advective or reversible part of transport) with mixing (i.e. the irreversible part of transport). In this context, “adaptive” means that not only do the positions of the APs within the grid evolve during advective transport along Lagrangian trajectories, but also the number of APs dynamically adjusts based on the flow. This adjustment, always applied between pure trajectory steps in CLaMS, enhances the spatial resolution of the grid when strong flow deformations occur, while simultaneously reorganizing the grid to keep the total number of APs within certain limits. In the following, we revisit the key concepts of the CLaMS adaptive grid, with a focus on its role in representing convection – a critical component of transport between the Earth's surface and the middle and upper troposphere.

2.1 Advection

To describe the 3D position r of an AP in CLaMS, longitude ℓ and latitude φ are used as the horizontal coordinates, while the hybrid potential temperature ζ represents the vertical coordinate (Mahowald et al., 2002). The notation $\mathit{r}(\mathrm{\ell },\mathit{\phi },\mathit{\zeta },t)$ denotes the exact position of the AP in space and time. “Hybrid” refers to the combination of the dry potential temperature θ, which remains constant during adiabatic motion, and the orography-following coordinate $\mathit{\sigma }=p/p_{\mathrm{s}}$ (where p is pressure and p_s is surface pressure). This combination enables the model to resolve transport processes in both the stratosphere and troposphere, where orography plays a role. The transition between σ and θ is controlled by the parameter ${\mathit{\sigma }}_r=p/p_{\mathrm{s}}$, set to 0.3 in CLaMS, with pressure p at location r. The Lagrangian trajectory, which quantifies the advective motion of an AP in CLaMS, is determined by numerically integrating the equation of motion:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\dot{\mathit{r}}:={\displaystyle \frac{\mathrm{D}\mathit{r}}{\mathrm{D}t}}=\mathit{u}(\mathit{r},t),\mathit{r}\left(t_{\mathrm{0}}\right)={\mathit{r}}_{\mathrm{0}},\mathit{u}=(u,v,\dot{\mathit{\zeta }})\end{array}$$

Here, u and v are the horizontal velocities (in m s⁻¹), and $\dot{\mathit{\zeta }}$ is the hybrid vertical velocity (in K d⁻¹), composed of the diabatic heating rate $\dot{\mathit{\theta }}$ and $\dot{\mathit{\sigma }}$, the latter being a function of the kinematic vertical velocity in pressure coordinates ($\dot{p}$) and the surface pressure tendency ($\dot{p_{\mathrm{s}}}$) (for details, see Mahowald et al., 2002; Ploeger et al., 2010). In CLaMS, the Runge–Kutta 4th-order method is used for the numerical integration of Eq. (1). The 4D fields u, v, $\dot{\mathit{\theta }}$, $\dot{\mathit{\sigma }}$, and $\dot{p}$ are provided by the reanalysis, here by ERA5 (Hersbach et al., 2020).

The black solid line in Fig. 1 shows an example of the time evolution of a typical Lagrangian trajectory calculated from prescribed 3D winds provided by the reanalysis. On timescales of hours to days, atmospheric motion is predominantly isentropic (i.e. with θ∼const.), not only in the stratosphere but also in large parts of the free troposphere, far from convective cells. Cross-isentropic motions mainly occur over longer timescales and are primarily driven by radiative imbalances (Krüger et al., 2008; Schoeberl and Dessler, 2011; Legras and Bucci, 2020). The exception is convection, which uplifts air from the Planetary Boundary Layer (PBL) into the free atmosphere, where latent heat release during the condensation and freezing of water vapor plays a dominant role. The red line in Fig. 1 illustrates such a trajectory, featuring a strong diabatic and highly localized (net) updraft Δθ, followed by subsequent slow, large-scale diabatic subsidence, mainly due to longwave cooling of radiatively active water vapor transported from the main convective outflow (e.g. Poujol and Bony, 2024 and references therein). If these convective updrafts are resolved by the meteorological wind data, they appear primarily in the vertical velocity fields $\dot{\mathit{\zeta }}$. If not, i.e. if unresolved, the respective updrafts must be parameterized. As described in Konopka et al. (2019, 2022), CLaMS triggers unresolved convection using specific values of either the Convective Available Potential Energy (CAPE) or the moist Brunt–Väisälä frequency, N_m. The value of Δθ is derived from the total available latent heat of the water vapor (transforming from gas to ice) in the CLaMS AP within the PBL, as provided by the reanalysis.

Figure 1 Longitude (ℓ), latitude (φ), and hybrid potential temperature (ζ) describe the spatial positions of air parcels (APs) in CLaMS. The trajectory quantifies the advective motion of an AP in a time- and space-dependent velocity field $(u,v,\dot{\mathit{\zeta }})$, which is updated at a frequency of $\mathrm{1}/\mathrm{\Delta }t$, typically with Δt ranging from 1 to 6 h. This is also the frequency at which unresolved convection – i.e. convection not captured in the velocity fields – is incorporated. Two trajectory examples are shown: a quasi-isentropic trajectory (i.e. almost adiabatic, θ∼const.) (black), and a strongly diabatic trajectory affected by a convective updraft (Δθ) (red). The convective updraft lifts the AP from the Earth's surface to a point above the surface θ=const., after which it gradually descends back to the θ=const. level. The red dotted points represent the projection of this trajectory onto the θ=const. surface.

2.2 CLaMS grid

The core concept of CLaMS leverages the numerical diffusion that arises from interpolations between APs to parameterize atmospheric mixing, as implemented in the adaptive grid procedure described in previous studies (McKenna et al., 2002a; Konopka et al., 2004, 2012). The numerical horizontal and vertical diffusion, D_h and D_v, associated with interpolations between two neighboring APs, scale as $\sim r^{\mathrm{2}}/\mathrm{\Delta }t$ and $\mathrm{\Delta }{z}^{\mathrm{2}}/\mathrm{\Delta }t$, where r and Δz represent the geometric horizontal and vertical mean separation between the APs, and $\mathrm{1}/\mathrm{\Delta }t$ represents the interpolation frequency (see, e.g. Press et al., 2007, Sect. 20.2: Diffusive Initial Value Problems). Consequently, the ratio $D_{\mathrm{h}}/D_{\mathrm{v}}={\mathit{\alpha }}^{\mathrm{2}}$, where $\mathit{\alpha }=r/\mathrm{\Delta }z$, quantifies the atmospheric aspect ratio, reflecting the ratio of atmospheric variability scales (Haynes and Anglade, 1997). Since this relation aligns with observations (Balluch and Haynes, 1997), it suggests that simple interpolations and their associated numerical diffusion can effectively mimic atmospheric diffusion. Thus, the parameters r and Δz, which are generally time- and space-dependent, can be used to define the Lagrangian grid. The distances between neighboring APs – i.e. the grid itself – and the frequency of adaptive re-gridding are therefore critical for accurately capturing the diffusive properties of the entire atmosphere. CLaMS relies on forward Lagrangian transport and mixing among APs that cover the free atmosphere, from the top of the PBL up to the mesosphere, following the concept described by Konopka et al. (2012). In our standard configuration, almost 2 million APs are distributed from the surface to the mesosphere. The grid is defined by two key parameters: entropy density s and the dry Brunt–Väisälä frequency N_d, from which the grid parameters r (horizontal separation) and Δz (vertical separation) are derived. First, we ensure that each AP represents the same amount of entropy s, which guarantees a consistent representation of mixing processes across different atmospheric regions, as mixing is an irreversible process that increases entropy. In a materially closed system, the specific internal energy depends solely on the specific entropy s, so APs with the same entropy also share the same internal energy. Second, we assume that the aspect ratio α∼N_d. This assumption helps shape the adaptive grid to reflect the atmosphere's varying mixing properties due to stability differences, with a well-mixed tropopause and a stably stratified stratosphere.

As air mass density (ρ) decreases and dry potential temperature (θ) increases with altitude, the entropy density (s∼ρln θ) reaches a maximum – using an idealized 1D US standard atmosphere – near the tropical tropopause. Consequently, an even distribution of s across all APs results in the highest grid point density around this region (see Fig. 2). Given that N_d peaks in the lower stratosphere, and following our assumption that α∼N_d, the corresponding APs exhibit large horizontal extensions and reduced vertical separations, consistent with the observed large values of α in this region (Balluch and Haynes, 1997). In contrast, APs in the middle troposphere and upper stratosphere exhibit significantly smaller α values (cf. the ratios of the horizontal to vertical edges of the blue rectangles in Fig. 2). This relationship between N_d and α allows us to mimic the diffusive properties of the atmosphere, which is characterized by a well-mixed troposphere and strongly suppressed vertical mixing in the stratosphere, where static stability N_d is very high (see e.g. Grise et al., 2010). Using idealized 1D profiles of entropy s and Brunt–Väisälä frequency N_d from the US standard atmosphere, the CLaMS grid, defined by r(ζ) and Δz(ζ) (or rather Δζ(ζ)), is generated with r=r₀ at ζ=380 K (tropical tropopause), based on ζ values provided by the US standard atmosphere with σ_r=0.3 (Konopka et al., 2012). By varying r₀, different model resolutions can be chosen, with the standard configuration set by r₀=100 km (see Appendix B).

Figure 2 CLaMS grid extending from the Earth's surface to the mesosphere (grey points). The density of air parcels is highest around the tropical tropopause. Four idealized APs are shown: three in the free atmosphere (blue) and one in the Planetary Boundary Layer (PBL, red). While the APs in the free atmosphere are generated following the 1D entropy-static stability concept described in Konopka et al. (2012), the generation of APs in the PBL is decoupled from this procedure to ensure sufficiently high resolution, in both time and space, to resolve small-scale processes such as convective updrafts not represented in the meteorological velocity fields or various sources of tracers.

While this grid design, based on simple thermodynamics, provides a good approximation for the free atmosphere, it becomes less applicable near the Earth's surface, where transport is strongly influenced by boundary effects such as friction and convection. In this region, interactions with the Earth's surface and PBL dynamics introduce additional small-scale sources of tracers and convective updrafts from small areas (∼1 km² and below), which are critical for accurate modeling of atmospheric transport. Examples of such small-scale sources include localized emissions from urban centers, industrial facilities, and biomass burning events, as well as biogenic emissions from vegetation. To account for these boundary effects, additional information must be integrated into the model. In information theory, Shannon entropy is sometimes used to quantify these additional entropy fluxes (Shannon, 1948). In our simplified approach, we decouple the APs in the PBL from those in the free atmosphere and significantly increase the resolution within the PBL by reducing the volume of the respective APs. This enhanced resolution more accurately captures the finer-scale processes characteristic of near-surface atmospheric dynamics.

After initializing the grid as described above, the advective trajectories transport the APs to their new positions. If necessary, the chemistry module updates the chemical composition of each AP. The advective time step, Δt, is typically set to 24 h, though shorter values between 6–24 h may also be used (note that the internal time step in the Runge–Kutta scheme for trajectory calculation is much smaller). Once the advection step is completed, the CLaMS adaptive grid procedure is applied to parametrize mixing and rearrange the grid.

2.3 Adaptive grid procedure

This procedure consists of two parts: the insertion of new APs when the nearest neighbors (NNs) of the considered AP move too far apart during the time step Δt (deformation-induced regridding) and the merging of too-close NNs to maintain the total number of APs within a prescribed range. The key numerical challenge in both parts is the computationally expensive calculation of the NNs, which, in the current CLaMS version, is effectively performed in 2D using Delaunay triangulation (see Appendix A). Since isentropic mixing typically occurs in relatively thin, quasi-2D atmospheric layers, this approach is physically well-justified. The adaptive grid procedure is thus performed layer by layer, dividing the full atmosphere into nearly isentropic layers i, with the altitude-dependent thickness Δζ_i defined by the entropy s and the static stability N_d profiles introduced in the previous subsection.

As shown in Fig. 3, deformation-induced regridding, which parameterizes isentropic mixing, considers only those APs whose final vertical positions (i.e. after the advection step) are within the layer Δζ_i and whose initial vertical positions are within a small but prescribed range Δζ, ensuring an isentropic approximation. An AP affected by convective uplift from the boundary layer during the advection time step, as exemplified by the red AP, is excluded from this part of the adaptation procedure because it significantly violates the isentropic approximation.

Figure 3 Layerwise adaptive grid procedure in CLaMS with unresolved convection (red AP). After the advection step, all APs are grouped into horizontal layers with thickness Δζ_i. For example, this applies to AP P and the other seven APs in the right ellipse (for better visibility, the layer Δζ_i is shifted downwards, with dashed lines marking the respective vertical positions of all considered APs). These seven APs are considered because they are the (old) nearest neighbors (NNs) of P before the advection step (left ellipse), i.e. their horizontal distance to P was ∼r, where r represents the mean separation of APs in layer Δζ_i. Since the vertical positions of these APs before the advection step lie within a small but prescribed range Δζ (typically Δζ≈Δζ_i), their respective advection is considered quasi-isentropic. The deformation-induced mixing (the first part of grid adaptation) inserts new APs if the distance between NNs becomes larger than r₊ (as is the case for segments PA, PB, and PC, marked by green crosses). In the second part of grid adaptation, APs are merged if their separation from their (current) NNs is smaller than r₋. Only during this part can APs that have moved from outside of Δζ into Δζ_i, such as the red AP denoting convection, undergo interpolation due to merging.

Consider an AP and its position before being advected by a Lagrangian trajectory. We calculate its NNs using the horizontal coordinates prior to advection, referred to as “old” coordinates, to distinguish them from the “current” coordinates after the advection step. Deformation-induced regridding means that new APs are inserted between P and its NNs (calculated with the old coordinates) if their separations r after the advection increase beyond a prescribed value r₊ (in our example, this applies to segments PA, PB, and PC), with relevant quantities interpolated at these new positions from their respective neighbors (for details on interpolation, see Appendix C). Note that r before the advection represents the mean horizontal distance in the considered layer Δζ_i. We refer to this part of grid adaptation as deformation-induced regridding, as the drifting apart of neighboring APs is mainly caused by vertical shear and horizontal strain, typically understood as precursors to mixing in stably stratified flow.

Since deformation-induced regridding increases the total number of APs (in our example, by 3), the merging step of the adaptive grid procedure is applied to balance the insertion of new APs and maintain the total number within a prescribed range. For this, a renewed 2D triangulation is performed using the current positions of the APs, and air parcels are merged if their distance is smaller than a prescribed threshold r₋. In this process, the red point introduced by convective uplift may also be merged with other APs. Note that the advection and convection time steps are decoupled in this approach: the advection time step is typically 24 h, while the convection time step corresponds to the maximum frequency of the available reanalysis data (every 6 h here, though 1 h data are also available). The convection time step is also typically the interval when the lower boundary, defining species values in the PBL, is updated with new data.

The choice of parameters r_± is somewhat arbitrary. While r₊ depends on Δt and should quantify the strength of critical deformations that trigger mixing, the approach where $r_{+}=r\exp \left(\mathit{\lambda }\mathrm{\Delta }t\right)$ (with r as the horizontal mean separation in layer i and λ the Lyapunov exponent) has proven effective for triggering AP insertion. The choice of r₋ is more pragmatic, motivated by the need to confine the total number of APs and prevent unlimited growth in the number generated during insertion. We use $r_{-}=\mathit{ϵ}r\exp (-\mathit{\lambda }\mathrm{\Delta }t)$, with ϵ ranging between 0.3–1, for Δt=6 and 24 h, respectively. This empirical choice of ϵ ensures that the total number of APs varies by less than 40 %. For typical values of the parameters Δt, λ, and ϵ, see Table 1.

Table 1 Key differences among the CLaMS versions. The Planetary Boundary Layer (PBL) is characterized by its thickness (Δζ) and the mean distance between air parcels (APs), denoted by r. The adaptive grid procedure is defined by its frequency (1/Δt), the Lyapunov exponent (λ), and the parameter ϵ, which collectively determine the criteria for the insertion and merging of APs. The parameterization of convective uplifts, representing unresolved convection, is triggered by either the moist Brunt–Väisälä frequency (N_m) or the Convective Available Potential Energy (CAPE), with critical values set to 0 for N_m in CLaMS-2.0 and 800 J kg⁻¹ for CAPE in CLaMS-3.0. In CLaMS-2.0, the adaptive and convective time steps are synchronized, whereas in CLaMS-3.0, they are decoupled. Furthermore, CLaMS-3.0 uses weighted interpolations during the insertion and merging of APs (see Appendix C), unlike the simple linear interpolations employed in earlier versions.

2.4 A Few Remarks

2.4.1 Layerwise adaptive grid and hybrid potential temperature

Ideally, the layers in the adaptive grid procedure should be parallel to the surfaces of constant potential temperature θ, ensuring “nearly” isentropic interpolations within the layers. When using the hybrid potential temperature ζ with σ_r=0.3, this condition is generally satisfied above about 300 hPa, away from strong orography. In regions influenced by orography, and up to approximately 200 hPa, ζ surfaces intersect the isentropic levels. The use of a hybrid ζ coordinate is particularly advantageous in this context, as it allows the model to better capture adiabatic motions such as gravity waves, which are less well handled by purely kinematic schemes. However, in regions of strong orographic forcing, this problem is still not completely solved. Increasing the parameter σ_r would make the model more diabatic, as the respective vertical velocities would then be calculated from the diabatic heating rates $\dot{\mathit{\theta }}$ rather than from the kinematic vertical velocity $\dot{p}$ (see Appendix A in Tao et al., 2018). The optimal parameter choice still requires further investigation.

2.4.2 2D Delaunay triangulation

The layerwise adaptive grid procedure uses 2D Delaunay triangulation to calculate the NNs of all APs in the layer, both before and after the advection step. 2D triangulation projects all APs onto a sphere and, in this way, neglects vertical separation between the APs. Since potential NNs across the boundary between the layers are neglected, this approximation is slightly weakened by considering alternating layer thickness every second advective step Δt (Konopka et al., 2004). 2D Delaunay triangulation on a sphere is not only numerically more efficient but also conceptually simpler, as there are no boundaries on a sphere.

2.4.3 Mass conservation

As in Eulerian models, a known concentration or mixing ratio, together with the volume of the grid box, provides the mass, which can also be determined for the Lagrangian grid, although the Eulerian grid cell must be replaced by the Voronoi volume (see Appendix A). Thus, the mass of a species, such as a chemically passive tracer, can be calculated both locally and for the entire atmospheric domain. It is important to note that interpolations within the CLaMS adaptive grid do not inherently conserve mass (see Appendix C). A few points are worth mentioning here. Regarding global mass conservation within the model, the total mass of each chemically passive species in the resolved atmospheric domain should balance with the fluxes into and out of this domain. This is trivially true in stationary (i.e. time-independent) atmospheric states, where all concentrations and fluxes are constant, and consequently, the system's mass remains unchanged. Thus, in cases where chemical production is roughly balanced by chemical destruction, and the atmospheric composition is solely driven by seasonality, mass conservation is inherently preserved in CLaMS as the model remains close to a stationary state.

Alternatively, Lagrangian grids can be formulated with a constant number of APs, each containing the same mass. While such grids maintain a constant total mass by design, they face significant challenges. First, the number of APs decreases with altitude as atmospheric mass declines exponentially with height. Consequently, fewer APs are available in the upper atmosphere, potentially insufficient to accurately capture relevant physical processes. Second, if the 3D velocity fields do not precisely satisfy the continuity equation, the initially uniform distribution of APs can become distorted, leading to local violations of mass conservation. For example, in such a model, all APs could migrate into one hemisphere, compromising the model's physical accuracy. The CLaMS adaptive grid addresses these issues by inserting new APs to fill potential gaps or reducing the number of APs in overly clustered regions. This approach ensures that the model maintains high-quality transport even in the presence of these challenges.

3 Adaptive Grid for Different CLaMS Versions

To visualize the performance of the CLaMS adaptive grid, we compare the standard configuration of CLaMS-3.0 with CLaMS-1.0 (Pommrich et al., 2014) and CLaMS-2.0 (Konopka et al., 2022). The key differences between these versions are summarized in Table 1. Three separate CLaMS runs are initialized on 1 January 1979 and span 44 years, reaching the end of 2023. The initial grid of APs is generated as described in the previous section, using the 1D entropy-static stability concept (Konopka et al., 2012). The nominal horizontal resolution is set to r₀=100 km, with an aspect ratio of α=250 at ζ=380 K. The thickness and horizontal mean distance between the APs in the planetary boundary layer (PBL) vary slightly between the versions, as outlined in Table 1, with Δζ ranging from 100 to 140 K. The choice of Δζ=140 K corresponds to a geometric thickness varying between 1.6–2.2 km.

During the simulation, after each advective time step Δt, the adaptive grid procedure adjusts the resolution of the initial grid based on flow deformations. Typically, the total number of APs grows by 30 %–40 % – from approximately 1.3, 1.4, and 1.8 million for CLaMS-1.0, 2.0, and 3.0, respectively – before stabilizing throughout the integration. Figure 4 illustrates the time evolution of globally integrated profiles for the total number of APs in CLaMS-3.0 (left panel) and the differences in such profiles for all three CLaMS versions, as exemplified on 31 August 2017 (right panel). The gray line in Fig. 4 represents an idealized number density distribution, assuming that all APs have the same mass. The CLaMS distributions generally follow this mass-scaling between 380–450 K but intentionally deviate from it above 500 K, where more APs are used to resolve the atmospheric structure in this region.

Figure 4 Vertical profiles of the globally integrated total number of APs, plotted as step-wise functions with step heights corresponding to the layer thicknesses. Left: Time evolution for the CLaMS-3.0 configuration. Right: Differences between the CLaMS versions on 31 August 2017. The gray line represents a distribution scaled by atmospheric density, assuming that each AP has the same mass.

After an initial spin-up period of approximately two months, the total number of APs in CLaMS-3.0 stabilizes at around 2.4 million, with little change thereafter (left panel). On 31 August 2017, a representative snapshot day, three distinct peaks in AP density are observed around 600, 340 K, and between 140–260 K. These peaks are present across all three versions of CLaMS (right panel), with CLaMS-3.0 extending to the Earth's surface in contrast to the previous CLaMS versions. Figure 5 further illustrates these features, showing zonal means of AP density (left) and mixing intensity (right), using an area-preserving latitude grid horizontally and the initial ζ-grid vertically. The three maxima in AP density correspond to the locations of jets and the main convective outflow in the tropics. Specifically, they highlight increased AP density around the polar vortex over Antarctica and the subtropical jets, as well as in the tropical main convective outflow. The mixing intensity reflects the percentage of APs influenced by regridding and numerical diffusion, with higher values near the jets indicating more active mixing. Additionally, a distinct feature in the tropics around 500 K is driven by the Quasi-Biennial Oscillation (QBO), where strong wind shear between easterlies above and weaker westerlies below causes flow deformations and drives mixing. Strong convective activity below the tropical tropopause moves APs from the PBL into the Tropical Tropopause Layer (TTL), enhancing AP density and triggering regridding in the adaptive grid procedure.

Figure 5 Adaptive grid properties for CLaMS-1.0, 2.0, and 3.0, shown for 31 August 2017 as an example. Zonal means of AP number density (n, left) and mixing intensity (Mix. Int., right) are presented, calculated per grid box defined by an area-preserving latitude grid and the ζ-grid generated at initialization. The mixing intensity quantifies the percentage of APs affected by the adaptive grid procedure. The vertical positions of the layers, total number of levels, and total number of APs are indicated on the right side of each panel. Solid and dashed lines represent westerlies and easterlies, respectively, while the thick white line represents the tropopause.

Similar features are observed near the tropospheric jets between 70–250 K. In CLaMS-3.0, higher AP density in the PBL further amplifies this effect. An asymmetry in stratospheric AP distributions between the Southern Hemisphere (SH) and Northern Hemisphere (NH) is evident, particularly in the CLaMS-1.0 and 3.0 configurations, reflecting the dynamic contrast between the vigorous winds of the winter SH and the calmer winds of the summer NH, sometimes described as “solid body rotation”. In contrast, this asymmetry is less pronounced in the CLaMS-2.0 configuration, which appears to be more affected by numerical noise, with a weaker representation of the QBO effect. This discrepancy arises from the higher frequency of grid adaptation in CLaMS-2.0 (6 vs. 24 h in CLaMS-3.0), making the detection of deformations more sensitive to numerical noise during the re-gridding. Despite the improved representation of tropospheric transport in CLaMS-2.0 compared to CLaMS-1.0 (Konopka et al., 2022), we exclude CLaMS-2.0 from further discussions and focus henceforth on comparing CLaMS-3.0 with CLaMS-1.0.

4 Model Evaluation with CO Satellite and In Situ Data

The most significant updates in CLaMS-3.0 pertain to the representation of upward transport from the PBL to the UTLS region. Carbon monoxide (CO) serves as an excellent tracer for validating the model's troposphere-to-stratosphere transport due to its relatively short lifetime of around 3 months, which enables the tracing of distinct signatures of convective transport (Park et al., 2009; von Hobe et al., 2021; Pan et al., 2024). Moreover, CO benefits from extensive satellite-based observations (Hegglin et al., 2021) and numerous local in situ measurements, which can be utilized both for defining the CLaMS lower boundary and for validating the transport and chemistry implemented in the model.

In particular, satellite observations from the Measurement of Pollution in the Troposphere (MOPITT) and the Atmospheric Infrared Sounder (AIRS) provide suitable lower boundary conditions for the model, capturing both the spatial and temporal variability of CO with high accuracy. Using back- and forward-trajectory calculations, we grid the MOPITT (2001–2012, version 4) and AIRS (2013–2020, version 6) data at a latitude-longitude-time resolution of 2, 6°, and 5 d, respectively, at ζ=200 K (corresponding to ∼500 hPa) (Pommrich et al., 2014). These particular time periods were chosen to maintain continuity with previous studies and datasets used in earlier phases of the model's development, ensuring consistency across analyses. This approach establishes the CLaMS lower boundary by overwriting all CLaMS APs with ζ values below 200 K after each advective time step, Δt.

4.1 Comparison of CLaMS Model to MLS and ACE-FTS Satellite Observations

To evaluate CLaMS-3.0, we use Microwave Limb Sounder (MLS) version 5.1 (Livesey et al., 2020) and Atmospheric Chemistry Experiment Fourier Transform Spectrometer (ACE-FTS) version 5.2 (Sheese et al., 2015; Boone et al., 2023; Sheese and Walker, 2023) satellite observations spanning the period 2005–2020. MLS, a limb-sounding instrument, provides extensive global coverage with dense daily observations but has a coarser vertical resolution of approximately 5–5.5 km for CO (Santee et al., 2017). In contrast, ACE-FTS, based on solar occultation measurements, offers higher vertical resolution (1–3 km for CO; Hegglin et al., 2009) but with much sparser spatial and temporal sampling, constrained to specific sunrise and sunset conditions along its orbit.

We analyze partial columns of CO (XCO), calculated in mg m⁻², for both satellites and CLaMS across the pressure range from 146.8 to 68.1 hPa, as discussed in Wright et al. (2025). This metric approximately captures the entire UTLS region in both tropical and extratropical zones, providing a convenient integral quantity to study the zonally resolved seasonality of CO. While the high vertical resolution of ACE-FTS allows for direct comparison with CLaMS, applying the MLS averaging kernels to the higher-vertical-resolution model data would be necessary for direct comparisons with MLS profiles. However, since the CLaMS partial columns calculated with and without applying the MLS averaging kernel differ by less than 5 % (see Appendix D), we restrict our comparisons to simple partial columns derived from both MLS and CLaMS. The zonally resolved annual mean climatologies for 2005–2020 for both satellites and CLaMS simulations, along with their zonal means, are shown in Fig. 6.

Figure 6 Annual mean climatology of CO partial columns (XCO) between 147–68 hPa, derived from MLS and ACE-FTS data for the period 2005–2020 (first two panels on the left), and from CLaMS-1.0, CLaMS-3.0, and CLaMS-3.0 with enhanced CO at the lower boundary (CLaMS-3.0/sens) (right). The respective zonal means are shown in the bottom left panel.

Figure 7 Seasonal anomalies (DJF, MAM, JJA, and SON) of the CO partial column, ΔXCO, derived from the 2005–2020 climatology for MLS (left), CLaMS-1.0 (middle), and CLaMS-3.0 (right). The anomalies are calculated relative to the zonal mean for the respective season.

Table 2 Seasonal Spatial Correlation Coefficients between Observed CO Anomalies (MLS/ACE-FTS) and Model Simulations (CLaMS-1.0 and CLaMS-3.0) as shown in Fig. 7.

Sensitivity studies for 2006, a year when these differences are most pronounced, show that incorporating a parameterization for unresolved convection is critical for reproducing such spatial anomalies (Table 3, second row). However, the choice of parameterization type (CAPE-based or the moist Brunt–Väisälä frequency N_m) has a smaller impact (not shown). The increased density of APs (40 vs. 110 km as in CLaMS-1.0) plays a much weaker role in representing anomalies (Table 3, third row) but is more crucial in enhancing the absolute values of CO columns (not shown).

Table 3 Sensitivity analysis for 2006: Spatial correlation coefficients between MLS and different CLaMS-3.0 configurations. The first row represents the reference case, where all sensitivities are switched on (1 = enabled, 0 = disabled). The other rows show the impact on correlation when each sensitivity is disabled.

First, as ACE-FTS and MLS climatologies act as references, it should be noted that they differ both in absolute values and spatial patterns, with MLS values being approximately 10 % higher in the tropics and more than 100 % higher poleward of 50°. As discussed in Livesey et al. (2020), a substantial high bias (factor of ∼2) in version 2 (v2) MLS CO, as revealed through comparisons with various in situ CO datasets, has been essentially eliminated in version 5 (v5). However, the quoted systematic uncertainties in v5 MLS CO remain unchanged from those in version 4 (v4), meaning the 2σ values provided for the UTLS in the study by Santee et al. (2017) are still valid. These uncertainties, which are rather large ($\pm \mathrm{30}\mathrm{ppbv}+\pm \mathrm{30}\mathrm{\%}$), may be sufficient to encompass the magnitude of the biases identified in the comparisons with ACE. On the other hand, zonal asymmetries compare well between the two climatologies, with one notable exception: the region over northwest Australia, where high CO values are visible in ACE-FTS but not in the MLS climatology. Here, sampling biases in ACE-FTS due to its non-uniform coverage cannot be ruled out (Kloss et al., 2019). Such differences make model validation somewhat challenging.

Our first robust finding is that the observed CO partial column values are significantly underestimated by CLaMS-1.0, while CLaMS-3.0 shows notable improvements. However, the model values of XCO remain smaller than those from the satellite instruments. As can be deduced from the comparison of zonal means (Fig. 6, bottom left), CLaMS-3.0 underestimates XCO from ACE-FTS by around 15 % in the tropics and only by a few percent elsewhere. The low bias relative to MLS is larger, especially in the extratropics, similar to the differences between ACE-FTS and MLS.

As vertical transport from the PBL to the UTLS might still not be fully captured in CLaMS-3.0, another potential limitation is the relatively coarse spatial resolution (2 and 6° in latitude and longitude) and temporal resolution (every 5 d) of the lower boundary derived from MOPITT and AIRS. These observations are smoothed by their averaging kernels, potentially leading to an underestimation of CO values. To assess this effect, we performed a sensitivity study (CLaMS-3.0 with enhanced CO) by scaling the lower boundary with factors of 1.5 and 1.25 for the MOPITT and AIRS datasets, respectively, reflecting the generally smaller CO values in MOPITT. Although somewhat arbitrary, these scaling factors were chosen to adjust the lower boundary to upper-range values typically observed in the PBL (see also Pommrich et al., 2014; Vogel et al., 2016). In the tropics, this adjustment reversed the low bias in CLaMS XCO to a slight high bias relative to MLS (see Fig. 6, bottom left panel). However, a significant low bias still persists in the extratropics relative to MLS.

Before introducing in situ data to investigate this discrepancy further, we compare seasonal anomalies between CLaMS-3.0 and MLS, which offers denser spatial coverage compared to ACE-FTS. The results are shown in Fig. 7, with panels from top to bottom representing DJF, MAM, JJA, and SON.

To assess the agreement between observed CO anomalies and those simulated by CLaMS-1.0 and CLaMS-3.0, we computed spatial correlation coefficients for each season using data from both MLS and ACE-FTS. The results are presented in Table 2. Notably, CLaMS-3.0 consistently shows higher correlation coefficients across all seasons compared to CLaMS-1.0, indicating a more accurate representation of the spatial distribution of CO anomalies. In particular, positive XCO anomalies over Africa during DJF and MAM, as well as over South America during DJF and SON, are more pronounced in CLaMS-3.0. Additionally, the enhanced influence of the Asian summer monsoon during JJA is better captured in CLaMS-3.0 compared to CLaMS-1.0. The superior spatial coverage of MLS data compared to ACE-FTS likely contributes to the higher correlation coefficients observed with CLaMS-3.0. This discrepancy is particularly evident during SON over northwest Australia, where elevated CO values are detected by ACE-FTS but are absent in both MLS climatology and our CLaMS-3.0 simulations.

Interestingly, replacing CO with an E90-tracer-based diagnostic (Prather et al., 2011; Abalos et al., 2017), which mimics CO behavior through a 90 d decay, and substituting the satellite-based lower boundary with a uniform E90 tracer value of 150 ppbv globally, also shows some potential to reproduce the anomalies (Table 3, last row). This highlights the relative importance of accurately representing vertical transport over the specific details of the lower boundary conditions. Moreover, the relatively high correlation with the E90 tracer further motivates numerical experiments with such an idealized tracer, as they are simpler in design and do not require a lower boundary condition. This simplification makes it easier to optimize transport schemes, for example, in climate models.

In Appendix E, we complement our validation of the UTLS model representation by analyzing ozone partial columns from both CLaMS-3.0 simulations and MLS observations. As ozone is predominantly a stratospheric tracer, unlike carbon monoxide (CO), this comparison primarily assesses the model's performance in the stratospheric part of the UTLS. CLaMS-3.0 reproduces the climatological seasonality and longitudinal structure of MLS-derived ozone remarkably well, with some improvement over CLaMS-1.0. Notably, the model now captures more pronounced ozone minima over the Maritime Continent, reflecting improved convective lofting of ozone-poor air – consistent with the CO-based diagnostics in Sect. 4. This confirms the continued ability of the model to accurately represent stratospheric transport processes (see, e.g. Konopka et al., 2010; Yan et al., 2018).

Figure 8 Observed in situ CO profiles from various campaigns in the Northern Hemisphere (gray points represent all local flights of the respective aircraft) and their mean values (black thick points) calculated over bins defined by the MLS pressure levels, with the respective 2σ standard deviations. The profiles are plotted as a function of log-pressure altitude, using H=7 km and a reference pressure of 1000 hPa. The corresponding mean profiles from CLaMS simulations are also shown (see legend). The MLS and ACE-FTS profiles represent mean values and their standard deviations over the months and regions covered by the respective campaigns. The thick black lines on the right y axes indicate the partial column range of XCO, spanning 147–68 hPa, calculated for the MLS and ACE-FTS partial columns.

Table 4 Summary of campaigns, including the geographic base, months covered, latitude-longitude ranges used in the MLS and ACE-FTS climatologies for comparison with in situ data (Fig. 8), and CO instruments utilized: COLD/COLD2 (Carbon monoxide Laser Detector) (Viciani et al., 2008, 2018), UMAQS (University of Mainz Airborne QCL-based Spectrometer) (Krasauskas et al., 2021), and Aerodyne-CO (Aerodyne Quantum Cascade Laser Absorption Spectrometer) (McManus et al., 2010; ACC, 2023).

4.2 Comparison with In Situ Data

To further investigate the discrepancies in XCO absolute values, especially the high bias of MLS relative to ACE-FTS and CLaMS-3.0 poleward of 50°, we incorporate in situ airborne observations (Fig. 8) from several campaigns. The RECONCILE (Reconciliation of Essential Process Parameters for an Enhanced Prediction of Arctic Stratospheric Ozone Loss and Climate Effects) campaign was carried out in 2010 onboard the M-55 Geophysica aircraft, based in Kiruna, Sweden (von Hobe et al., 2013). The StratoClim (Stratospheric Climate Links with Emphasis on the Upper Troposphere and Lower Stratosphere) campaign was conducted in 2017, also onboard the M-55 Geophysica aircraft, based in Kathmandu, Nepal (von Hobe et al., 2021). The WISE (Wave-driven Isentropic Exchange) campaign, also in 2017, operated from Shannon, Ireland, using the HALO (High Altitude and Long Range Research Aircraft) (Lauther et al., 2021). The ACCLIP (Asian Summer Monsoon Chemical and Climate Impact Project) campaign took place in 2022, based in Osan, South Korea, using the GV (Gulfstream V) and WB-57 aircraft (Pan et al., 2024). The PHILEAS (Probing High Latitude Export of Air from the Asian Summer Monsoon) campaign was conducted in August and September 2023, with most flights operating from Anchorage, Alaska, onboard HALO.

The campaigns, the geographic positions of their bases, the months covered by the observations, and the longitude-latitude ranges assumed for the comparison with the MLS and ACE-FTS climatologies are listed in Table 4. The CO datasets collected during these campaigns are shown in Fig. 8, using log-pressure altitude as the vertical coordinate, and are compared with the respective profiles derived from the MLS and ACE-FTS climatologies, as well as CLaMS simulations. The comparison reveals that, in agreement with our discussion in the previous subsection, MLS profiles consistently show higher values than ACE-FTS profiles. For instance, MLS XCO values exceed those from ACE-FTS by more than 100 % in the polar region (RECONCILE), with the in situ mean profile exhibiting better agreement with ACE-FTS observations than with MLS observations. In other campaigns – except for WISE – particularly above 10 km, the mean in situ observations generally agree better with MLS profiles.

The largest deviation between the in situ data and satellite observations was observed during the ACCLIP campaign, likely due to the sampling of highly polluted Asian monsoon air and the strong La Niña influence discussed in Pan et al. (2024), which makes the climatologies unrepresentative for this campaign. CLaMS-3.0 simulations generally fall within the spread of in situ observations (within the 1σ range around the mean profiles) but tend to slightly underestimate the observed values. The best agreement is achieved when using the enhanced lower boundary configuration, except for the StratoClim comparison, as relatively clean air was sampled within the Asian monsoon anticyclone during that campaign. Only CO observations onboard the Geophysica (RECONCILE and StratoClim) and WB-57 aircraft (ACCLIP) cover the full range of XCO considered in this study. However, only the RECONCILE observations clearly show better agreement with ACE-FTS than with MLS, further confirming the high bias of MLS CO measurements northward of 50° N.

5 Conclusions

This study presents an updated and improved version of the Chemical Lagrangian Model of the Stratosphere (CLaMS-3.0/MESSy), featuring key advancements in the representation of transport processes from the Planetary Boundary Layer (PBL) to the Upper Troposphere and Lower Stratosphere (UTLS) region. By enhancing spatial resolution in the PBL and decoupling the frequency of the adaptive grid procedure from the convection parameterization, CLaMS-3.0 more accurately captures the complex interplay between small-scale convective uplift and large-scale flow dynamics. These changes allow for better representation of critical transport pathways and mixing processes, especially for tracers like carbon monoxide (CO).

Using CO as a tracer, we evaluate CLaMS-3.0 against Microwave Limb Sounder (MLS) and Atmospheric Chemistry Experiment Fourier Transform Spectrometer (ACE-FTS) satellite observations, as well as airborne in situ measurements from multiple campaigns. Including a parameterization for unresolved convection significantly improves the representation of CO anomalies in the UTLS region, defined here as the partial columns between 147–68 hPa. Seasonal patterns over Africa, South America, and Asia are now captured more effectively, leading to better alignment with satellite climatologies. Our results confirm that MLS-based partial columns consistently show higher CO values compared to ACE-FTS, with differences largest in the extratropics poleward of 50°. CLaMS simulations align more closely with ACE-FTS than with MLS in that region, and in situ observations in the boreal polar region generally support ACE-FTS rather than MLS. Conversely, in situ measurements agree better with MLS than with ACE-FTS in air affected by the Asian summer monsoon above 10 km.

Sensitivity tests with enhanced boundary-layer CO demonstrate the importance of accurately representing lower boundary conditions for CO and the potential impact of unresolved sources in this region. Incorporating such improvements would further reduce discrepancies with in situ and satellite observations, particularly in convectively active regions. This underscores the need for better spatial and temporal resolution of boundary-layer emissions to enhance model accuracy globally. The improvements in CLaMS-3.0 are not merely incremental; they mark a significant step forward in the model's ability to resolve the intricate transport and mixing processes in the UTLS region. For instance, the enhanced resolution in the PBL and the decoupling of grid adaptation from convection parameterization enable more realistic simulations of convective outflows by applying unresolved convective updrafts every 6 h and mixing every 24 h. In future work, we plan to increase the frequency of convection parameterization to 1 h to further enhance model performance.

In summary, CLaMS-3.0 offers a robust and refined tool for studying UTLS transport processes. The advancements in grid resolution and parameterization of unresolved convection make the model better suited for future studies exploring seasonal and spatial variability in the UTLS region, the impacts of boundary-layer emissions, and the representation of convective transport in a changing climate.

Appendix A: Delaunay triangulation

Handling irregular grids, especially determining the nearest neighbors (NNs) in such grids, requires specific mathematical approaches. The most intuitive definition of the NNs for a given point in an irregular grid embedded in a k-dimensional space can be achieved using the concept of the Voronoi volume (or Voronoi area in 2D). The Voronoi volume is defined as the region where every point is closer to the considered point than to any other grid point. In Fig. A1, the light pink-colored region marks the Voronoi area of an arbitrary point in an example 2D grid.

The borders of the Voronoi regions are called Voronoi edges (dashed lines) and represent the set of points equidistant from two neighboring grid points. In this way, both the volume and the NNs of the considered grid point are defined. Note that Voronoi edges also serve as symmetry axes for pairs of neighboring grid points. Furthermore, every grid point lies within exactly one associated Voronoi region. Two grid points that are nearest neighbors are connected by a line called a Delaunay edge which forms a network of triangles in 2D, known as the Delaunay triangulation, or in 3D, the Delaunay tetrahedration.

Determining the Delaunay triangulation directly (i.e. without using the concept of Voronoi regions) for an arbitrary irregular grid in k dimensions is a challenging numerical problem. One approach, the convex hull method (Barber et al., 1996), with computational costs scaling as Nlog N (N being the total number of grid points), is based on the idea that finding the Delaunay triangulation for a set of points in k-dimensional Euclidean space can be converted into the problem of finding the convex hull for a set of points in (k+1)-dimensional space (Preparata and Shamos, 1985). The convex hull of a set of points is the smallest convex set containing this set. For a planar set of points, the convex hull can be visualized by imagining an elastic band stretched around the outermost points; when released, it will assume the shape of the required convex hull. (Note that due to the appropriate projection to a higher dimension, in this example from 1D to 2D, no points will be “left” in the apparent interior of our example 2D set of points).

It can be shown that the concept of layerwise calculation of the NNs in a 3D domain containing N grid points, as used in the CLaMS adaptive grid procedure, is less time-consuming than a full 3D triangulation, which scales as N² in the worst case (de Berg et al., 2008). Let us consider m horizontal layers fully filling the 3D domain, with each layer consisting of n² grid points. Thus, the aspect ratio $\mathit{\alpha }=r/\mathrm{\Delta }z$, where r and Δz are the mean horizontal and vertical separation between the grid points, scales as $\mathit{\alpha }=n/m$, and the total number of grid points can be rewritten as $N=m{n}^{\mathrm{2}}=n^{\mathrm{3}}/\mathit{\alpha }$. Consequently, the ratio of the scaling times between the full 3D triangulation, t_3D, and the layer concept, $t_{{\text{2D}}^{\prime }}$, can be written as:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle \frac{t_{\text{3D}}}{t_{{\text{2D}}^{\prime }}}}\sim {\displaystyle \frac{N^{\mathrm{2}}}{m{n}^{\mathrm{2}}\log n^{\mathrm{2}}}}\sim {\displaystyle \frac{n^{\mathrm{6}}}{n^{\mathrm{3}}\log n^{\mathrm{2}}}}>{\displaystyle \frac{n^{\mathrm{6}}}{n^{\mathrm{5}}}}=n\sim N^{\frac{\mathrm{1}}{\mathrm{3}}}\end{array}$$

where we used the estimate log x<x, and n was expressed in terms of N. Thus, as N increases, the 3D triangulation becomes increasingly time-consuming. Furthermore, the 3D triangulation of grid points at the boundaries of a typical atmospheric 3D domain, which extends from the Earth's surface to the upper edge of the atmosphere, is much more complicated compared to the 2D triangulation on closed spherical layers without such boundaries.

Figure A1 Voronoi area (light pink) confined by the Voronoi edges (dashed lines) of an arbitrary grid point within an irregular 2D grid, also defining the volume associated with this point and its nearest neighbors (NNs). Delaunay triangulation (bold lines) is a direct method for finding the NNs.

Appendix B: Entropy-static stability grid

Figure B1 2D vs. 1D idealized forms of the entropy density s (top) and the dry Brunt–Väisälä frequency $N_{\mathrm{d}}^{\mathrm{2}}$ (BVF, bottom). The 2D functions represent the ERA5 2005–2020 all-year zonal means (left, black lines are isotachs of zonal winds), while the 1D profiles (right) are derived from the meridionally averaged 2D distributions. For comparison, the 1D profiles based on the US standard atmosphere, as discussed in Konopka et al. (2012), are also shown.

Following the procedure described in Konopka et al. (2012), the altitude z dependent grid parameters r (horizontal separation) and Δz (vertical separation) are given as:

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\displaystyle }r={\left({\displaystyle \frac{kN{S}_{\text{ap}}}{s}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\mathrm{\Delta }z={\left({\displaystyle \frac{\mathrm{1}}{kN}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}{\left({\displaystyle \frac{S_{\text{ap}}}{s}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\end{array}$$

with $S_{\text{ap}}=r_{\mathrm{0}}^{\mathrm{2}}\mathrm{\Delta }{z}_{\mathrm{0}}{s}_{\text{max}}$ denoting the entropy per air parcel with horizontal and vertical extensions r₀ and Δz₀. The 1D profiles of the entropy density s and the dry Brunt–Väisälä frequency N_d are derived from the US standard atmosphere and are shown in the top and bottom right panels of Fig. B1. The index “0” refers to the altitude of the entropy maximum, s_max, where $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{0}}=r_{\mathrm{0}}/\mathrm{\Delta }{z}_{\mathrm{0}}=\mathrm{250}$ is the expected aspect ratio in the lowermost stratosphere around θ=380 K (Balluch and Haynes, 1997). The respective value of $N^{\mathrm{2}}(\mathit{\theta }=\mathrm{380}\mathrm{K})=N_{\mathrm{0}}^{\mathrm{2}}=\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{s}}^{-\mathrm{2}}$ allows us to determine the parameter $k={\mathit{\alpha }}_{\mathrm{0}}/N_{\mathrm{0}}$. By varying r₀, i.e. the model horizontal resolution at θ=380 K (roughly corresponding to the position of the tropical tropopause), different global model resolutions can be chosen, with r>r₀ and Δz>Δz₀ above and below the level where the 1D entropy s has its maximum (see Fig. 4 in Konopka et al., 2012).

However, the question arises as to how this idealized 1D approach can be justified by more realistic spatial and temporal distributions of s and N. Although a full 3D approach would be desirable, the climatological all-year zonal means of these two functions, shown in the left part of Fig. B1, may offer some insights into an optimal grid structure.

Following the definition of entropy density, $s=c_p\mathit{\rho }\ln (\mathit{\theta }/{\mathit{\theta }}_{\mathrm{0}})$ (in J (K m³)⁻¹), with specific heat c_p (J (K kg)⁻¹) and air mass density ρ (kg m⁻³), the ERA5-based, all-year 2D zonal mean of s, as a function of latitude and altitude, is shown in the top left panel of Fig. B1. The freely chosen reference potential temperature θ₀ is set to 300 K, as the lowest potential temperature surface that does not cross the Earth's surface. Thus, s values below this level are negative and cannot be used to calculate the grid parameters r and Δz from Eq. (B1). By zonally averaging this 2D distribution, the corresponding 1D profile is shown in the top right panel of Fig. B1, together with the 1D profile derived from the US standard atmosphere (Konopka et al., 2012) (note that for the 1D case θ₀ is set to the potential temperature at the Earth's surface).

As the 2D s distribution, being a product of the upwardly decreasing air mass density ρ and the increasing potential temperature θ, shows distinct maxima in the polar regions around 15 km and in the tropics around 5 km, the corresponding 1D distribution has a broad maximum around 13–14 km. Consequently, the highest density of APs will be generated in regions different from those suggested by the 2D distribution. Additionally, the 2D s distribution implies that grid generation should be decoupled not only in the PBL but also in the entire region where θ<θ₀. Here, θ≥θ₀ is considered the free atmosphere, and θ<θ₀ as an “extended” PBL. Note that the isentropes below θ₀ cross the Earth's surface and, assuming isentropic motion, are inherently connected with the PBL. Similar discrepancies in the vertical positioning of maxima can be observed between the 2D and 1D $N_{\mathrm{d}}^{\mathrm{2}}$ profiles, as shown in the bottom panels of Fig. B1.

Nonetheless, the 1D profiles of s and $N_{\mathrm{d}}^{\mathrm{2}}$ derived from the US standard atmosphere approximate fairly well those derived from the ERA5 climatology. Additionally, while the seasonality of the 2D distributions may play a role (not shown), it is less pronounced for s than for $N_{\mathrm{d}}^{\mathrm{2}}$. Using the 2D distribution to generate the grid would be more physically accurate and desirable, as it could better capture spatial and temporal variability. However, implementing such an approach faces unresolved challenges, such as handling layers that are no longer disjunct, which complicates the grid structure. As a result, the 1D approach remains the practical choice in the current implementation of CLaMS.

Appendix C: Adaptive grid interpolations

Consider two APs with their respective (Voronoi) volumes V₁ and V₂. The number of molecules N_i of a trace gas within each AP can be expressed as:

$$\begin{array}{}\text{(C1)}& {\displaystyle }{\displaystyle }{N}_i={\mathit{\mu }}_i{n}_i{V}_i,i=\mathrm{1},\mathrm{2}\end{array}$$

where μ_i and n_i denote the mixing ratio of the trace gas and the air number density, respectively, for each AP. Assume that these air parcels undergo a mixing event in CLaMS, requiring the mixing ratio of the resulting mixed AP to be derived from the mixing ratios of the original APs.

For the resulting volume V and air density n of the mixed AP, we assume:

$$\begin{array}{}\text{(C2)}& {\displaystyle }{\displaystyle }V=V_{\mathrm{1}}+V_{\mathrm{2}}\text{and }nV=n_{\mathrm{1}}{V}_{\mathrm{1}}+n_{\mathrm{2}}{V}_{\mathrm{2}}.\end{array}$$

Using algebra, the mixing ratio μ of a trace gas in the mixed AP can be written as:

$$\begin{array}{}\text{(C3)}& {\displaystyle }{\displaystyle }\mathit{\mu }={\displaystyle \frac{N_{\mathrm{1}}+N_{\mathrm{2}}}{nV}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{1}}{n}_{\mathrm{1}}+\mathit{\eta }{\mathit{\mu }}_{\mathrm{2}}{n}_{\mathrm{2}}}{n_{\mathrm{1}}+\mathit{\eta }{n}_{\mathrm{2}}}}\end{array}$$

where $\mathit{\eta }=V_{\mathrm{2}}/V_{\mathrm{1}}$. Since both APs lie within the same CLaMS layer, their volumes are approximately equal by construction at initialization. Over time, the adaptive grid procedure may lead to slight inhomogeneities in the AP distribution within a layer due to atmospheric deformation. However, we neglect this effect here and assume η=1 (to approximate η, 2D Voronoi areas, as discussed in Appendix A, can be applied). Note that determining the exact 3D volumes of Lagrangian APs is challenging, making other assumptions for η difficult to implement. After some simplification, Eq. (C3) becomes:

$$\begin{array}{}\text{(C4)}& {\displaystyle }{\displaystyle }\mathit{\mu }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\mathit{\mu }}_{\mathrm{1}}(\mathrm{1}+\mathit{\chi })+{\mathit{\mu }}_{\mathrm{2}}(\mathrm{1}-\mathit{\chi })\right]\end{array}$$

with $\mathit{\chi }=\mathit{\delta }/n$, $n=\mathrm{0.5}\cdot (n_{\mathrm{1}}+n_{\mathrm{2}})$, and $\mathit{\delta }=\mathrm{0.5}\cdot (n_{\mathrm{1}}-n_{\mathrm{2}})$. Using the ideal gas law $n\sim p/T$, we derive:

$$\begin{array}{}\text{(C5)}& {\displaystyle }{\displaystyle }\mathit{\chi }={\displaystyle \frac{\mathit{\delta }}{n}}={\displaystyle \frac{\frac{p_{\mathrm{1}}}{T_{\mathrm{1}}}-\frac{p_{\mathrm{2}}}{T_{\mathrm{2}}}}{\frac{p_{\mathrm{1}}}{T_{\mathrm{1}}}+\frac{p_{\mathrm{2}}}{T_{\mathrm{2}}}}}.\end{array}$$

Appendix D: Impact of the MLS Averaging Kernel

Figures 6 and 7 present comparisons of CLaMS and MLS CO partial columns (XCO) calculated over the UTLS range from 147 to 68 hPa for both datasets. While the MLS data are provided at three pressure levels – 147, 100, and 68 hPa – with a vertical resolution of approximately 3 km, their true vertical resolution (quantified by the full width at half maximum of the rows of the averaging kernel matrix) is about 5–5.5 km at these altitudes (Livesey et al., 2020). In contrast, the CLaMS data have a much higher vertical resolution, on the order of a few hundred meters. To enable an appropriate comparison between MLS data and CLaMS output, the averaging kernel provided by MLS can be applied to smooth the model data vertically, mimicking the coarser vertical resolution of the satellite. Figure D1 quantifies this effect on both the partial columns and their anomalies.

Figure D1 Annual mean climatologies of the CO partial column (XCO, in mg m⁻²) for the period 2005–2020, spanning from 147 to 68 hPa (left) and their respective zonal anomalies (ΔXCO, right). Results are shown from CLaMS interpolated to six pressure levels (top), to three MLS pressure levels (middle), and with the MLS averaging kernel applied to these three levels (bottom).

In the top row of Fig. D1, both quantities are calculated from CLaMS output interpolated onto six pressure levels: 146.8, 121.2, 110.0, 100.0, 90.0, 82.5, 75.0, and 68.0 hPa, which approximately match the vertical resolution of the Lagrangian grid. In the second row, both XCO and its anomaly are calculated on the native MLS CO levels in the UTLS, i.e. 146.8, 100.0, and 68.0 hPa. Additionally, the MLS CO averaging kernel (which is provided on MLS levels) is applied, with the results displayed in the bottom panels of Fig. D1. While the reduction from six levels to three levels increases XCO values by no more than 2.2 %, the smoothing with the averaging kernel reduces the column by a maximum of 5 %. The effect on the anomalies is even smaller, by an order of magnitude. Consequently, we do not apply the MLS averaging kernel to the CLaMS results in this paper, as its impact is negligible.

Appendix E: Partial Columns of Ozone

Figure E1 Seasonal (DJF, MAM, JJA, and SON) values of the O₃ partial column, XO₃, derived from the 2005–2020 climatology for MLS (left), CLaMS-3.0 (middle), and CLaMS-1.0 (right; shown as the percentage difference relative to CLaMS-3.0).

To further validate the model's representation of the UTLS region, we complement the CO analysis by evaluating partial columns of ozone (XO₃), defined as ${\mathrm{XO}}_{\mathrm{3}}={\int }_{p_{\text{bot}}}^{p_{\text{top}}}({\mathrm{O}}_{\mathrm{3}}/g)\mathrm{d}p$, where O₃ is the ozone mass mixing ratio, g is the gravitational acceleration, and the integration is carried out between p_bot=147 hPa and p_top=68 hPa. This pressure range aligns with native Microwave Limb Sounder (MLS) levels and thus facilitates a direct comparison with satellite data.

XO₃ values are derived from CLaMS-3.0/1.0 simulations and compared with MLS observations. As with CO, we find that the integral values are only weakly affected by whether the MLS averaging kernel is applied to the model output (difference below 10 %, not shown).

Figure E2 Same as Fig. E1, but for O₃ partial column seasonal anomalies, ΔXO₃, calculated for CLaMS-3.0 relative to the zonal mean for each respective season.

While Fig. E1 shows the seasonality of the absolute values, the anomalies are shown in Fig. E2. The correlation coefficients between XO₃ absolute values from CLaMS-3.0 and MLS range from 0.82 (SON) to 0.99 (JJA), with a value of 0.95 for the multi-seasonal average (ALL), and are even higher for the anomalies. While the overall improvement compared to CLaMS-1.0 is minor in terms of correlation, the difference patterns, CLaMS-1.0 minus CLaMS-3.0 (last column in Fig. E1), show a clearer signal in the tropics. In particular, over the Maritime Continent, CLaMS-3.0 captures stronger lofting of ozone-poor marine boundary layer air into the UTLS, consistent with the enhanced role of convection and in line with the CO-based results discussed in Sect. 4.

These findings confirm that the model developments presented here primarily influence the tropospheric part of the UTLS, where CO – rather than ozone – is a more suitable tracer. Since the model's lower ozone boundary (approximating the PBL at 1.6–2.2 km) is set to zero and tropospheric ozone chemistry is neglected (Pommrich et al., 2014), we refrain from overinterpreting the absolute XO₃ values. Nevertheless, as a stratospheric tracer, ozone remains well represented, with CLaMS-3.0 retaining the key features already captured by CLaMS-1.0 (Konopka et al., 2010; Yan et al., 2018), while extending model improvements into the lower UTLS.