© The Authors 2025

1 Introduction

The formation and growth of planets starts around a young star during the protoplanetary disc phase. Although the most direct observations of protoplanets so far (Keppler et al. 2018; Haffert et al. 2019; Mesa et al. 2019; Currie et al. 2022; Hammond et al. 2023) are of (super-)Jupiter-mass protoplanets orbiting at tens of astronomical units (au), the inner few au are also of special interest, being the location where inward pebble drift leads to a higher midplane dust-to-gas ratio (Weidenschilling 1977), which may be conducive to planetesimal formation (Drążkowska et al. 2016). The inner disc is also the site of the water snow line – where dust at the disc midplane becomes hot enough for H₂O to thermally desorb and enter the gas phase –, which may assist the concentration of these building blocks and subsequent planet formation just outside the snow line (Stevenson & Lunine 1988; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Schoonenberg et al. 2018). More massive planets may also potentially form in the inner disc, or otherwise form at larger radii, where substructures often attributed to planets are commonly observed (Andrews 2020; Bae et al. 2023), and migrate inwards (Lodato et al. 2019). As protoplanets will inherit their compositions from the disc material, information on the chemistry of the inner disc region sheds light on the initial chemical abundances of planets that form in, or migrate to, the inner disc.

These inner regions are warm enough to excite line emission at the mid-infrared (MIR) wavelengths, which has allowed the volatile molecular content of this planet-forming material to be studied with several generations of high-resolution spectrographs and space-based telescopes (Pontoppidan et al. 2014). JWST is now providing an unprecedented view of this gas-phase chemistry (van Dishoeck et al. 2023; Kamp et al. 2023; Henning et al. 2024), and its intriguing diversity across the disc population. For example, the multitudinous abundant hydrocarbons and weak or absent water seen towards very low-mass stars (Tabone et al. 2023; Arabhavi et al. 2024) confirms the carbon-rich nature of their discs previously inferred from abundant C₂H₂ in Spitzer spectra (Pascucci et al. 2009, 2013). On the other hand, in T Tauri stars with masses of ≳0.3 M_⊙, the main molecules tend to be H₂O and CO₂. In this case, diversity is seen between ‘CO₂-dominated spectra’, such as those of GW Lup (Grant et al. 2023) and CX Tau (Vlasblom et al. 2025) in which the Q-branch of CO₂ is particularly prominent (Pontoppidan et al. 2010) – and in which the ¹³CO₂ isotopologue may also be detected –, and spectra where the forest of water lines is the most prominent, such as those of CI Tau, GK Tau, HP Tau, IQ Tau (Banzatti et al. 2023), Sz 98 (Gasman et al. 2023), SY Cha (Schwarz et al. 2024), and DR Tau (Temmink et al. 2024b,a).

To understand the origin of this diversity, we must identify the factors that determine the inner disc chemistry: for example, we must decipher whether this chemistry reflects the chemical equilibrium, given the temperature, stellar radiation field, dust distribution, and cosmic ray flux (e.g. Eistrup et al. 2016, 2018), or a non-equilibrium maintained by significant transport of volatiles, either in the gas phase following the accretion flow or as ices on inwardly drifting dust grains. Moreover, the extent to which the observed chemistry relies on which of the layers of the disc are probed by the measurements remains uncertain.

To see how chemistry and transport compete, Booth & Ilee (2019) compared a static disc model to models with purely viscous evolution and then also to models incorporating dust radial drift (see also Cevallos Soto et al. 2022). For example, for CO₂ or NH₃, Booth & Ilee (2019) find that the timescale for gasphase viscous transport inwards from their snow lines is shorter than the timescale on which the bulk of the material at the midplane becomes transformed by gas-phase chemical reactions in the inner disc. This allows transport to entirely counteract any decrease in their abundances in the inner disc, which is agreement with the results of Bosman et al. (2018), who find that highly enhanced destruction rates would be needed to significantly depress inner disc CO₂ abundances. On the other hand, Booth & Ilee (2019) also find that for a more volatile molecule such as CH₄, which is transformed (by conversion into larger hydrocarbons) on a shorter timescale relative to its viscous transport time (given the low typical viscosities inferred at tens of au; e.g. Rosotti 2023), viscous transport is unable to erase the signature of chemical reactions. Radial drift further enhances the speed of the inward transport of the volatiles. This in turn further increases the degree to which major volatiles such as H₂O and CO₂ are enhanced, and negligible differences are seen in their profiles between models with and without chemical reactions. The strong CH₄ abundance decrease in the inner disc is now prevented by replacement by drifting ices, despite an efficient conversion to hydrocarbons when reactions are included.

Consequently, disc evolution models have been used to explore the evolution of the C/H, O/H, and C/O ratios in the inner disc due to the successive enhancement of H₂O, CO₂, and then finally hydrocarbons (such as CH₄) and CO by drifting pebbles (e.g. Booth et al. 2017). This happens because the ices of these molecules have increasingly large binding energies and therefore have snow lines at increasingly large radii. This means that these molecules are deposited by desorption from drifting pebbles successively further from the star, and then take longer to be advected into the inner disc by the gas. However, the absolute timing of these phases may vary, for example as a function of stellar mass, as lower-mass stars have correspondingly cooler discs, closer-in snow lines, and shorter timescales for transport (Mah et al. 2023).

However, most studies which investigate these phenomena neglect the fact that the aforementioned substructures are frequently observed and are typically understood to represent dust traps that would interrupt the supply of volatiles and could solve several issues in planet formation and disc evolution modelling (Birnstiel 2024). While focused on photo-evaporation, Lienert et al. (2024) also include a single example of a planet-carved gap in a model of disc chemical evolution by transport, but do not systematically study the effect of the gap’s properties. A more systematic exploration – albeit only in relation to H₂O – was conducted by Kalyaan et al. (2021, 2023), who demonstrated the impact of dust traps at different radii on the H₂O vapour mass in the inner disc. A key finding was that, the closer the trap to the star, the greater the fraction of the H₂O ice reservoir it could isolate from reaching the inner disc and hence the more the enhancement of the H₂O vapour mass due to drift is suppressed. Moreover, Kalyaan et al. (2021) argue that drift simultaneously shrinks dust discs and delivers water to the inner disc, and may thus explain the anti-correlation between the luminosity of midinfrared (MIR) water lines and the radius containing 90% of the continuum flux at millimetre wavelengths (Banzatti et al. 2020); however Kalyaan et al. (2021) do not directly forward model such quantities. Finally, Mah et al. (2024) also explored the role of the trap depth and formation time on the inner disc H₂O abundance and also the C/O ratio.

In order to understand the connection between the disc structure at large radii and the diversity in the MIR molecular spectra of T Tauri discs, we built on these previous works by modelling the delivery of both CO₂ and H₂O to the inner regions of discs with dust traps. We focus on these two molecules as they are (a) very widely detected in spectra from JWST’s MidInfrared Instrument Medium-Resolution Spectroscopy mode (MIRI-MRS) spectra (unlike e.g. NH₃) and are (b) relatively insensitive to chemical processing throughout the bulk of the disc (unlike hydrocarbons such as CH₄ or C₂H₂). We conducted this investigation using a 1D disc evolution code (Booth et al. 2017), and in Sect. 2 we summarise its key features and the modifications we made to the chemistry module. To make the connection to observations more direct, we then produced simple forward models of the MIR spectra and retrieved representative column densities, temperatures, and emitting areas using 0D slab models assuming local thermodynamic equilibrium (LTE) as described in Sect. 3. Section 4 presents the results of this exercise in terms of the evolution of the vapour masses and retrieved column densities. We compare the predicted ranges to observed discs, and explore any dependence on gap location at different stages of disc evolution. In Sect. 5 we discuss other processes or structural considerations that may impact our findings and consider other possible tests of our results. Finally, we summarise our conclusions in Sect. 6.

2 Disc evolution model

We used the 1D disc evolution code described in Booth et al. (2017); Booth & Clarke (2018); Sellek et al. (2020). The code implements the two-population dust model of Birnstiel et al. (2012) which approximates the dust distribution as a population of small grains that follow the gas and a population of large grains that dynamically decouple from the gas. The two populations have fixed mass fractions f_small and f_large = 1 − f_small respectively. These mass fractions are calibrated against fragmentation-coagulation simulations that account for two main barriers to growth which set the maximum grain size: a) the fragmentation of grains that collide with sufficiently high relative velocity and b) the removal of dust grains by radial drift faster than they can grow. These suggest f_small = 0.25 when fragmentation limits growth, or f_small = 0.03 when radial drift limits growth. The fragmentation limit usually applies when dust is abundant – so mainly while discs are still very young or where dust is concentrated – while more mature discs tend to be drift limited.

The initial disc surface density for initial mass M_D,0 and radius R_C,0 was given as a Lynden-Bell & Pringle (1974) profile: $${\text{Σ}}_{t=0}(R)=\frac{M_{\text{D},0}}{2\pi R_{\text{C},0}}\frac{1}{R}\exp \left(-\frac{R}{R_{\text{C},0}}\right).$$(1)

The temperature profile is assumed to follow a power law: $$T=T_0{\left(\frac{R}{R_0}\right)}^{-1/2},$$(2)

normalised such that the aspect ratio (scale height) at R₀ = 1 au is h₀ = 0.021 (H₀ = 0.021 au). This corresponds to a disc temperature at 1 au of T₀ = 118 K as assumed by Kalyaan et al. (2021) based on the median stellar luminosity in Banzatti et al. (2020), and also gives aspect ratios at larger radii in line with the M_* = 1 M_⊙ radiative transfer models of Sinclair et al. (2020).

In the following subsections, we describe more details of our approach to modelling the viscosity, dust, and traps. For ease of comparison, the fiducial parameters, which are summarised in Table 1, mostly follow Kalyaan et al. (2021, 2023).

2.1 Key features of gas and dust transport

The gas evolves due to the action of an effective viscosity, for which we assumed a Shakura & Sunyaev (1973) model v_eff = αC_SH with constant α. This results in a viscous velocity, which may be expressed in terms of the Keplerian velocity υ_K as $$\left|v_{\text{visc}}\right|=\frac{3}{2}\alpha h^2{v}_{\text{K}}.$$(3)

The dust interacts with the gas through drag; this interaction may be characterised by the Stokes number $$\text{St}=\frac{\pi }{2}\frac{{\rho }_{\text{s}}a}{{\text{Σ}}_G},$$(4)

where ρ_s is the dust material density and Σ_G is the gas surface density. This interaction results in the dust radial velocity having two components (e.g. Dipierro et al. 2018): one resulting from advection with the radial motion of the gas (see Eq. (3)), and another from radial drift (Whipple 1973) which is $$\left|v_{\text{rd}}\right|=\mathrm{St}{h}^2\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}R}{v}_{\text{K}}.$$(5)

As both 3/2 and the logarithmic pressure gradient $\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}R}$ are typically of order unity, the parameter St/α determines the coupling of the dust to the gas. For St/α ≳ 1, the radial drift term dominates, while for St/α ≲ 1, the viscous advection term dominates.

As discussed, in the early phase of disc evolution when the dust is abundant, or likewise in traps when there is no radial drift, the Stokes number of the largest grains is limited by fragmentation (Birnstiel et al. 2012). As collisions are typically most important between grains of similar sizes, with similar bulk velocities, the largest relative velocities are usually the random velocities due to turbulence, which results in the following limit: $${\text{St}}_{\max }=\frac{u_{\text{f}}^2}{3{\alpha }_{\text{t}}{c}_{\text{S}}^2},$$(6)

where u_f is the fragmentation velocity and typically it is assumed that α_t = α. Thus, for some typical values, the corresponding coupling parameter of the large, drifting dust grains in the vicinity of the water snow line at 150 K becomes $$\frac{{\text{St}}_{\max }}{\alpha }~0.6{\left(\frac{u_{\text{f}}}{1\text{\hspace{0.17em}m}{\text{\hspace{0.17em}s}}^{-1}}\right)}^2{\left(\frac{\alpha }{{10}^{-3}}\right)}^{-2}{\left(\frac{T}{150}\right)}^{-1}.$$(7)

We note, however, that if the bulk turbulence (as parametrized by α_t) is low, other processes can dominate the relative velocities (Pinilla et al. 2021). For example, fragmentation can occur due to the differential drift velocity between two dust grains of different sizes (Birnstiel et al. 2012), in which case the limit becomes $${\text{St}}_{\max }=\frac{u_{\text{f}}{v}_{\text{K}}}{\gamma c_{\text{S}}^2(1-N)},$$(8)

where $\gamma =\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}R}$, and N = 0.5 is the relative Stokes number assumed for the two differentially drifting grains.

In this work, we considered three possible scenarios for the parameters appearing in Eqs. (3)–(8) which we explain below and summarise in Table 2 and Figs. 1 and 2. These explore different possible regimes for the turbulence and collision velocities at which fragmentation limits dust growth, two key uncertainties in current dust evolution modelling (Birnstiel 2024).

1. The first parameter set assumes that icy dust grains are stickier and less fragile than dry grains, that is u_f,dust = 1 m s⁻¹ (Blum & Wurm 2008), u_f,ice = 10 m s⁻¹ (Gundlach & Blum 2015) and α_t = α = 10⁻³. In this case, Eq. (7) predicts that dry dust grains in the inner disc that have sublimated all of their ices will be coupled to the gas, while icy dust grains at large radii will decouple. This will lead to a ‘traffic jam’ (Saito & Sirono 2011; Pinilla et al. 2016, see also left-hand panel of Fig. 2) where the velocity of the grains decreases inside the water snow line and thus they pile up.

2. Secondly, we considered u_f,dust = 1 m s⁻¹, u_f,ice = 1 m s⁻¹, and α = 10⁻⁴. This is motivated by studies suggesting lower fragmentation velocities for icy grains (Gundlach et al. 2018; Musiolik & Wurm 2019); by also decreasing α, we maintain the degree of coupling in the outer disc, allowing effective radial drift to continue, while also preventing the traffic jam at small radii. We note, however, that lowering α will also slow down the gas evolution.

3. Finally to avoid simultaneously changing the gas evolution, we also investigated a case where α_t ≠ α (c.f. Pinilla et al. 2021, Model 6): u_f,dust = 1 m s⁻¹, u_f,ice = 1 m s⁻¹, α = 10⁻³ and α_t = 10⁻⁵ which should still allow the dust to form decoupled, drifting, pebbles both inside and outside the water snow line. In this case the fragmentation limit is due to the relative drift velocities (Eq. (8)).

Fig. 1 Cartoon depicting the three scenarios explored in this work. The thickness of the yellow arrows connecting the disc and star represents the magnitude of α controlling the gas evolution and accretion. The turbulence strength is indicated by the thickness and size of the circular arrow in the bottom right. The speckled region represents the dust with the brown colour indicating dry dust and the light blue colour indicating icy dust and the finer-grained pattern indicating where grains are more fragile. The vertical size of this region is inflated for illustrative purposes in locations where is a traffic jam and grains pile up.

Fig. 2 Example surface density profiles of gas (solid, initial value only) and (small+large) dust (dashed, at different times) for each of the three scenarios, illustrating the changing dust surface density distribution due to radial drift. The dip in the gas surface density signifies the gap, in this case at Rgap = 30 au and after 0.1 Myr (red), dust has begun to accumulate in the trap at its outer edge. The traffic jam that occurs in Scenario 1 can be seen in the steep rise in the dust surface density after 0.1 Myr inside the water snow line at ≈0.4 au. The other small steps correspond to changes in the solid fraction at snow lines and the smaller bumps are the accumulation of ices at snow lines due to the cold finger effect.

2.2 Traps

2.2.1 Imposition and shape of traps

We introduced dust traps as a perturbation to the effective viscosity v_eff – which governs bulk gas transport – as follows: $$v_{\text{eff}}=f(R,t)\alpha c_{\text{S}}H,$$(9)

where f (R, t) represents the additional torques needed to create and sustain a perturbation in the gas, for example due to the action of a planet. We note that these additional torques may not be turbulent or diffusive in nature, and so the radial diffusivity is assumed to be unaffected by the perturbation and as such expressed by D = αc_SH/Sc (where Sc ≈ 1 for St ≪ 1, is the Schmidt number) (Stadler et al. 2022; Kalyaan et al. 2023; Mah et al. 2024).

For the spatial dependence of f (R, t) we used a simple Gaussian defined in terms of the peak-to-baseline contrast A, the centre R_gap and the width w_gap: $$f(R,t)=1+(A-1)𝑔(t)\exp \left(-\frac{{\left(R-R_{\text{gap }}\right)}^2}{2w_{\text{gap }}^2}\right).$$(10)

Figure 2 includes an example of the gas surface density profile that results for A = 10 and R_gap = 30 au.

Observationally, the most easily quantified parameter of any trap is its radial location in the disc, so we explored a range of gap locations: R_gap ∊ {1, 3, 7, 15, 30, 60, 100} au. This spans from 1 au – in order to explore the impact of traps that may be too close to the star to resolve with the Atacama Large Millime- ter/submillimeter Array (ALMA) – to 100 au – just outside our fiducial value of R_C = 70 au. We note that the 1 au case will lie between the H₂O and CO₂ snow lines for the temperature profile we have assumed.

In order to be hydrodynamically stable, the width of the trap should be at least a scale height w_gap ≥ H (R_gap) (Lin & Papaloizou 1993; Dullemond et al. 2018). For consistency we follow Kalyaan et al. (2021, 2023) and used $w_{\text{gap}}=\sqrt{2}H\left(R_{\text{gap}}\right)$ (noting that our definition of w_gap differs from theirs by a factor $\sqrt{2}$). To ensure that traps of this width are resolved with several grid elements requires small $\frac{\text{Δ}R}{w_{\text{gap}}}$. Since (in general) $\frac{\text{Δ}R}{w_{\text{gap}}}\le \frac{\text{Δ}R}{H\left(R_{\text{gap}}\right)}=h_0^{-1}{\left(\frac{R_{\text{gap}}}{R_0}\right)}^{-1/4}\frac{\text{Δ}R}{R_{\text{gap}}}$, and the flaring is small (i.e. R^1/4 is a fairly weak function of R, so ${\left(\frac{R_{\text{gap}}}{R_0}\right)}^{-1/4}\sim 1$), then a sufficiently small $\frac{\text{Δ}R}{R}<h_0$ should suffice at all radii. Hence a logarithmically spaced grid is the natural choice for this case, rather than the more computationally efficient (for solving the viscous diffusion equation) ΔR ∝ R^1/2 which provides insufficient resolution at small radii and was thus found in testing to make close-in traps artificially leaky.

Previous work by Stadler et al. (2022) suggests that A ≳ 3 is needed to effectively trap pebbles, though somewhat lower values will still perturb the dust density profile. We used a fiducial A = 10 in order ensure we are comfortably in this regime, without creating very steep gradients that slow down the computation. The gap depth is commonly related to the equivalent mass of planet that would open that gap (e.g. Kanagawa et al. 2017). However, some caution is needed about making such interpretations. Firstly, the Kanagawa et al. (2017) parametrization represents the viscous criterion, namely how quickly the planet can clear the gap relative to the viscous process that act to refill it; however, an alternative limit is the thermal criterion (Lin & Papaloizou 1993), whereby the Roche radius of the planet should exceed the disc gas scale height.

Moreover, in models that solve the full Smoluchowski (1916) equation for coagulation (e.g. DUSTPY; Stammler & Birnstiel 2022), dust traps are known to be leaky (Drążkowska et al. 2019; Stammler et al. 2023): the small grains produced by collisional fragmentation in the trap are dragged with the gas – or more importantly diffuse – across the gap. Although we do not model the production of small grains by fragmentation in detail, as we employ a fixed f_small , our traps will still have a population of small dust that follows the concentration of the large dust grains. Thus the two-population model will produce diffusive leaking through the gap with a flux proportional to the diffusivity and the concentration gradient of the small dust: $$F_{\text{leak }}=2\pi R_{\text{gap }}D\text{Σ}\frac{\partial {ϵ}_{\text{small }}}{\partial R}.$$(11)

However, by construction, the fixed f_small of the two-population model does not account for this diffusive flux in determining how much of the dust is ‘small’ and susceptible to leaking and therefore our models will underpredict the proportion of the dust that is vulnerable to leaking. This caveat is not prohibitive to us modelling the transport of volatiles to the inner disc – so long as we are mindful of the fact that our gaps are potentially only producing the dust fluxes expected from a somewhat deeper gap. Consequently we refrain from reading too much into how the depth affects our results, or trying to equate our gap depth to a certain planet mass.

Finally, we also considered a smooth disc with no gaps at all to provide a reference model. Moreover we note that for simplicity, we do not consider discs with multiple gaps; the inner disc molecular abundances will generally be controlled by the innermost trap to the star (Kalyaan et al. 2021) as it is the final barrier isolating the dust – the inner disc is largely unaffected by whether the dust is actually trapped there or in a trap further out.

2.2.2 Growth of traps

One open question is at what stage of protoplanetary disc evolution most gaps and traps form, which is closely tied to the question of when planet formation begins. In order to see whether the effects of these structures on inner disc chemistry could be used to explore these questions, we explored a range of formation timescales as well as our fiducial case where the traps are present for all time (when we simply have 𝑔(t) = 1).

If the traps are produced by planets, they can only emerge once the planets have grown sufficiently to perturb the disc (e.g. Lin & Papaloizou 1993). When the additional torques represented by f (R, t) have a planetary origin, they scale with the square of the planet mass 𝑔(t) ∝ M_P(t)² (Lin & Papaloizou 1979). Meanwhile, planetary mass increases due to its accretion from the disc, usually modelled as some mass-dependent power law defined by ${\dot{M}}_{\text{P}}\propto M_{\text{P}}^{b_{\text{grow}}}$. Typical values are in the range 0.5 < b_grow < 3 depending on whether gas, pebble, or planetesimal accretion dominates and on what lengthscales the planet can accrete (Drążkowska et al. 2023). Thus, we adopted the following solution for the time dependent term in the torque: $$𝑔(t)=\min \left[1,{\left(1+\left(1-b_{\text{grow }}\right)\frac{t-t_{\text{f}}}{t_{\text{grow, f}}}\right)}^{2/(1-b)}\right].$$(12)

The free parameters in this model are the growth scaling b_grow, the time at which the torque ceases to grow t_f and the final (i.e. at t = t_f) growth timescale, t_grow,f.

The most important aspect of the gaps is when the perturbation to the disc profile grows large enough to trap pebbles. In the case of a planet, this is then associated with the isolation of the pebble flux feeding its growth. The final stage of pebble accretion typically involves accretion from the whole height of the disc (2D regime) with the relative velocities dominated by Keplerian shear for larger-Stokes-number pebbles (Hill regime) (Drążkowska et al. 2023). In this case, the growth scales as b_grow = 2/3 and can achieve a final growth timescale t_grow,f ~ 10⁴ yr at 1 au. Thus for our trap growth models, we fixed these values, and focussed on varying the end time t_f in the range t_grow,f < t_f ≤ few × Myr resulting in values t_f ∈ {0.03,0.1,0.3,1.0} Myr. Substructures are highly common in Class II discs in the best-studied regions with ages ≳Myr, and are also sometimes seen during the embedded phase, for example in HL Tau (ALMA Partnership et al. 2015) and IRS 63 (Segura-Cox et al. 2020) with ages ≲0.5 Myr. Our range of formation times thus straddles the appearance of the earliest known substructures.

2.3 Chemical tracers

The code includes a number of molecules as tracers in both gaseous and solid (ice) phase (aside from the refractory grain components, which are only ever solid). The gaseous tracers are advected with the gas velocity, while the ices move with the dust grains. For simplicity, in this work only transport (advection and diffusion) and state changes (freeze-out and desorption) are allowed to change the tracer abundances; there are no chemical reactions between species (which we validate in Appendix A). The only sink of a species is if it crosses the inner radial boundary of the model, at which point it is considered to have been advected onto the star.

2.3.1 Elemental abundances

The elemental abundances of planet-forming material are somewhat uncertain depending on what they are assumed to follow.

For example, the solar value (most recently estimated at C/H = 2.9 × 10⁻⁴ Asplund et al. 2021) is often used but is generally found to be higher than measurements of nearby young stars or the interstellar medium (ISM) suggesting the Sun is somewhat carbon-enriched, perhaps as a result of its formation location or environment, or galactic chemical evolution. For example, measurements of B stars (Przybilla et al. 2008; Nieva & Przybilla 2012), which are short-lived and therefore representative of contemporary star-formation, show median values of C/H = 2.1 × 10⁻⁴. Similarly, ISM measurements lie even lower at C/H ≈ 1.4 × 10⁻⁴ (Cardelli et al. 1996; Savage & Sembach 1996; Sofia et al. 2004) though this is understood as being due to some C being locked up in the solid phase. Nevertheless, recombination lines in the Orion Nebula suggest a higher gas-phase abundance of C/H ≈ 2.3 × 10⁻⁴ (Simón-Díaz & Stasinska 2011), which is suggestive of the return of the refractory components to the gas phase (Nieva & Przybilla 2012) in regions undergoing star formation. As this is intermediate between the B-star and solar inferences, we adopted this as an intermediate estimate of the total (volatile+refractory) elemental carbon budget.

However, the Orion Nebula appears depleted by an order of magnitude in silicates relative to solar, comparable to the ISM. Therefore, we cannot make the same assumption that refractory silicon carriers have sublimated as for the refractory carbon material. In this case, as cosmic abundances are highly consistent with solar, we simply adopted Si/H = 3.2 × 10⁻⁵ (Nieva & Przybilla 2012; Asplund et al. 2021) for the total Si abundance, all of which is assumed to be refractory.

Finally, the oxygen abundance resulting from the sum of the Orion Nebula value 4.5 × 10⁻⁴ and that which is locked up in silicates given the assumed silicate abundance and a 3:1 ratio is ~0.9 × 10⁻⁴ is 5.3 × 10⁻⁴. This value is intermediate between the present solar value (4.90 × 10⁻⁴; Asplund et al. 2021) and both the protosolar value (5.62 × 10⁻⁴; Asplund et al. 2021) and the B-star estimate of cosmic abundances (5.75 × 10⁻⁴; Nieva & Przybilla 2012); as it is consistent with both sets of values within the uncertainties, we adopted a total (volatile+refractory) O/H = 5.3 × 10⁻⁴.

Our adopted abundances are summarised in Table 3.

2.3.2 Molecular abundances

Correctly assigning the total elemental budgets to their various carriers is a non-trivial job. Adding up ices and expected volatiles typically falls comfortably short (e.g. Boogert et al. 2015; McClure et al. 2023); however, when including refractory sources one has to be simultaneously mindful of species that are very hard to constrain without distinct spectral features (such as amorphous carbon), and the risk of double counting species as not all spectral features are unambiguously assigned.

We assumed that the difference between our adopted carbon abundance and the ISM value is representative of the total refractory carbon, which suggests that around 40% of the carbon is in refractory material of some kind. For the volatile C-carrying molecules, we base our abundances on ice observations towards protostellar envelopes and molecular clouds (Öberg et al. 2011; Boogert et al. 2015). This only totals around 30% of the expected volatile carbon budget; we assume the missing C resides in CO, and add this to its contribution from the ices. The resulting total CO abundance is close to the total estimated from gas and ices for low-mass protostars (van Dishoeck et al. 2021).

For oxygen, we adopted an assumption that 20% of the total oxygen resides in H₂O. This is a factor 2–3 lower than expected by models (Bergin et al. 2000; Walsh et al. 2015) which predict highly efficient and rapid formation of water by grainsurface chemistry. Moreover, it is up to a factor 2–3 higher than derived for ices in protostellar envelopes or molecular clouds (Pontoppidan et al. 2004; Boogert et al. 2015; McClure et al. 2023). Dust growth is one plausible contributor to the low values: ices on larger grains, which could account for around 50% of the water ice, can become invisible to infrared absorption (van Dishoeck et al. 2021). However, as there is no direct evidence for this water, for the sake of tracking the missing oxygen in the code, we assume that the remaining 32% of the oxygen resides in O₂. Even though this is somewhat larger than predicted by thermochemical disc models (Walsh et al. 2015), and in tension with observed upper limits of a few % (van Dishoeck et al. 2021), this does not affect our results as we do not model its chemistry or emission.

2.3.3 Molecular data

Our chemical model is thus based on the Case 2 scenario of Booth et al. (2017); however, we include some additional species, updated some of the initial abundances and binding energies, and used prefactors in the desorption rates derived from experiments rather than the approximation of Hasegawa & Herbst (1993). The updated values for each species are listed in Table 4.

For the most abundant ices central to this work such as H₂O and CO₂, we expect many monolayers to be present on dust grain and thus we take our binding energies from multilayer Temperature Programmed Desorption experiments (which all correspond to pure ices). Only once the ice depletes to the final monolayer(s) will the binding energies become affected by the substrate and begin to deviate; although this will have little effect for the majority of the ice budget. For consistency, the binding energies are used in conjunction with prefactors derived simultaneously which is important for correctly identifying the location of the snow line (Minissale et al. 2022). These data are all taken from works included in Tables 2, 3, and 5 of the review by Minissale et al. (2022) with full details given in the notes of Table 4.

Thus, an implicit assumption is that the ices are pure and successively layered, each with a single desorption temperature. However observations of protostellar envelopes find that several components are needed to fit ice absorption bands, implying that the molecules may exist in several ice phases (see review by Boogert et al. 2015). In particular, for CO₂, up to five components are needed (Pontoppidan et al. 2008; Brunken et al. 2024), of which the pure phase is never seen to constitute more than about 20% (and is more typically <10%), while the polar phase (mixed with H₂O) is overall the most abundant (Öberg et al. 2011; Brunken et al. 2024). If all CO₂ resided in the polar phase, it could become trapped below frozen material even above its pure-phase desorption temperature. However, such CO₂ does not only codesorb at the H₂O desorption temperature, but much of it diffuses to the surface and sublimates; the rest is released at a temperature in between the CO₂ and H₂O sublimation temperatures in a process of ‘volcano desorption’ taking place during the amorphous-to-crystalline H₂O phase transition (Collings et al. 2004; Edridge et al. 2013). Moreover, thermal processing can lead to clusters of pure CO₂ forming from the polar phase in a process known as segregation (Ehrenfreund et al. 1999; Öberg et al. 2009), shifting the balance towards pure ices on longer timescales. Thus, given the likely mix of ice phases, possible segregation from the polar phase, and the release of CO₂ below the H₂O sublimation temperature even in the polar phase, we nevertheless expect a separation between the locations where H₂O and CO₂ sublimate in the disc; this work assumes the most extreme scenario.

3 Synthetic spectrum model

3.1 The origin of MIR emission

For many molecules, rotational, vibrational, and ro-vibrational transitions produce line emission in the mid-infrared. For CO₂, the strongest feature – often the only one clearly detected – is the Q-branch at 15 µm corresponding to purely vibrational excitation of the bending mode. In more CO₂ -dominated spectra, it can be flanked by the R- and P-branches from ro-vibrational excitation of the bending modes and the ‘hot bands’ at 13.85 µm and 16.2 µm – the Q-branches for transitions between bending and stretching modes. H₂O shows purely rotational transitions at wavelengths ≳10 µm, and ro-vibrational transitions with higher excitation energies (thus tracing hotter gas) at shorter wavelengths. These lines have typical excitation temperatures ≳1000 K and so trace regions with temperatures of 100s K. Thus, the most strongly emitting regions would be in the vicinity of the disc’s optical photosphere, where the temperatures rise, and/or the inner ≲1 au. This can be seen in 2D thermochemical models (Woitke et al. 2018; Bosman et al. 2022a,b), where the majority of line emission is from a thin hot layer. This layer is exposed to UV and undergoes rapid photochemical processing, and is not necessarily representative of the bulk composition towards the midplane (Anderson et al. 2021) that we follow with our 1D models. However, vertical transport and mixing may still lead to a signature in the surface of any midplane enhancement, depending on how the gas is reprocessed.

As well as the simple case of excitation determining the dominant emitting layer, there are three further features that limit the depth to which emission can be seen. Firstly, the dust in the disc provides a continuum opacity, which hides material beyond the MIR photosphere. Moreover, if there is a drop in the abundance vertically, for example when crossing the snow surface, emission may be potentially visible down to that depth but no further due to a lack of molecules to emit. Finally, the emission of strongest transitions may become optically thick due to the line opacity after a relatively small column density of material and no longer grow with an increased molecular column density; in this case changes in dust continuum opacity have less effect (Bruderer et al. 2015). This cooler material between the hot surface and the photosphere or snow surface may, in principle, be traced with an appropriate choice of lines, by picking those with a lower excitation energy which may be penalised by the partition function at high temperatures, or those which have low Einstein coefficients and thus can avoid becoming optically thick in the hot surface layer (e.g. Gasman et al. 2023). However in practice, optically thin lines may be rare (Temmink et al. 2024a, who found four clean, non-blended lines out of a total of >100, all tracing the hottest gas).

Moving beyond individual lines, the most widely used tools for the interpretation of MIR spectra are 0D LTE slab models. These are characterised by a column density, temperature, and emitting radius (itself a parametrization of the emitting area which simply scales the spectrum to match the observed fluxes) (e.g. Carr & Najita 2011; Salyk et al. 2011). Although many works start by fitting a single slab model to the emission from each molecule, 2–3 components are typically needed to fully represent the H₂O rotational spectrum, fitted either simultaneously using multiple slabs (Temmink et al. 2024a), to different parts of the spectrum (e.g. Gasman et al. 2023), or to the residuals after another disc spectrum is subtracted as a template (Banzatti et al. 2023); this allows weaker emission components to be recovered.

Typically, the dominant rotational H₂O component – and the one which a single slab finds - is the hottest, with T ≈ 600–800 K. This must represent the hot disc surface described and indeed matches the temperatures predicted for this layer by Bosman et al. (2022a) when chemical heating is accounted for. The coldest component, present only in some discs, is found at ≲200 K and is typified by prominent emission lines at 23.817 and 23.896 (Banzatti et al. 2023). Given its proximity to the desorption temperature of water, it may represent sublimating ices and has thus been interpreted as a signal of drift. However, emission purely in the vicinity of the midplane snow line would not reproduce the large emitting radius inferred (Banzatti et al. 2023; Temmink et al. 2024a) when these lines are prominent. Thus, this component must trace the full 2D nature of the snow surface, something that cannot easily be captured in our 1D models which assume a vertically isothermal disc. Between these two components is an intermediate ‘warm’ component, typically found at 250–500 K. This is the component that our 1D models are likely to mostly closely match in the absence of the full 2D thermochemical structure needed to understand the other components (see also Appendix C.2). However in practice, all its clean emission lines (i.e. detectable above noise, without blending with other features) may still be considerably optically thick (τ ~ 350, Temmink et al. 2024a) and hence not trace the full column density to the continuum MIR photosphere.

3.2 Constructing the spectrum

We follow an approach similar to Anderson et al. (2021) to construct a forward model of the MIR spectrum of our molecules of interest – H₂O and CO₂ – by dividing the disc into an annulus at each modelled radius R_i and summing the individual spectrum of each annulus F_λ,i, weighted by the area of the annulus: $$4\pi d^2{F}_{\lambda }={\displaystyle \sum _{i}2}\pi R_i\text{Δ}{R}_i{F}_{\lambda ,i},$$(13)

where d = 140 pc is a typical distance to a nearby disc, and where we sum over all annuli with a non-negligible gaseous reservoir of the element (here taken as all grid cells where f_spec,gas /f_spec,ice > 10⁻⁵). We modelled the spectrum over the wavelength range of MIRI-MRS Channels 3B (13.34–15.57 µm) & 3C (15.41–17.98 µm) which cover the main CO₂ features and the warm water emission (see Appendix C.2 for longer wavelength H₂O emission).

The spectrum of each annulus was given by a 0D LTE slab model calculated using RADEXPY (Tabone et al. 2023). As discussed, these are characterised by a temperature T_i and a visible column density of material N_vis,i. The slab models naturally account for the line opacity, so to model the emission of the bulk material, we must provide properties corresponding to the MIR τ = 1 surface which is the approximate maximum depth to which we can trace. Consequently, we took the midplane temperature T_i = T(R_i) given by Eq. (2) assuming the deep layers are approximately vertically isothermal. Assuming that gas – including the emitting molecules of interest – is vertically well mixed, then the molecules follow a Gaussian distribution in the vertical direction, and the visible column density is given by $$N_{\text{vis},i}=\frac{N_{\text{tot},i}}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{e}^{-{\tau }_i(\widetilde{ɀ})}}{e}^{-{\widetilde{ɀ}}^2}\text{d}\widetilde{ɀ},$$(14)

where $\tilde{z}=z/\sqrt{2H_i},{\tau }_i(\tilde{z}=\infty )=0,{\tau }_i(\tilde{z}=-\infty )={\tau }_i$, and N_tot,i is determined directly from the surface density of each molecule and its molecular mass. The main source of opacity is the small dust, which is also likely to be well-mixed (while the large grains largely settle to the midplane where they cannot contribute), and at each radius is assumed to have an abundance given by f_small ϵ where f_small = 0.03−0.25 comes from the two-population model, and the total dust-to-gas ratio ϵ varies spatially and temporally following the dust evolution that we model. In the absence of radial drift, this would imply a dust-to-gas ratio in the upper layers of 3 × 10⁻⁴–2.5 × 10⁻³, but we emphasise this may be raised by the delivery of dust or lowered if it drifts.

Given that the small dust is well-mixed and roughly follows the same vertical profile as the gas, we can understand Eq. (14) with the approximation that: $$N_{\text{vis},i}\approx N_{\text{tot},i}\frac{1-e^{-{\tau }_i}}{{\tau }_i}.$$(15)

In the limit τ_i ≫ 1, the visible column density is simply the fraction of the column lying above the MIR τ = 1 surface N_vis,i ≈ N_tot,i /τ_i.

We note, however, that for full accuracy in our calculation of τ(ɀ), we did include both the large and small grains, and calculate their vertical profiles include settling according to Dubrulle et al. (1995). For the opacity of the dust we use the opacities from DALI (Bruderer 2013) which represent a graphite-silicate mix as per the ISM composition of Weingartner & Draine (2001).

The most major limitation is that we assume a vertically constant temperature. In reality, different molecules are seen to present quite different emitting temperatures, with colder temperatures thought to trace lower heights into the disc. In general, CO₂ is expected to trace a deeper layer (i.e. higher column) than H₂O (Woitke et al. 2018), and ¹³CO₂ even deeper (and colder) due to excitation and line optical depth effects. Moreover, in equilibrium thermochemical models, H₂O is formed at a higher temperature than CO₂ which favours a high H₂O/CO₂ at high temperatures. A vertical temperature gradient would overall limit the column of material that dominates the emission, especially for the H₂O that traces the hotter material, and may alter the CO₂/H₂O ratio. In general, however, the height of optical photosphere which sets the depth of the hotter gas is also controlled by the dust column density, so the general relationship of Eq. (15) still holds, just with a higher optical depth.

For CO₂, we accounted for mutual shielding of lines (Tabone et al. 2023) – known as line overlap – but this is not required for the purely rotational lines of H₂O (Temmink et al. 2024a).

3.3 Interpreting the spectrum

To interpret our predicted spectra, we compared them to a grid of slab models covering T = 25−1500, K in steps of 25 K and ɪog₁₀(N/cm⁻²) = 16−25.92 in steps of 0.16 (following Gasman et al. 2023). As the flux of a slab model simply scales with the emitting area $A=\pi R_{\text{eff}}^2$, for each pair (T, N), the corresponding effective emitting radius R_eff is that which minimises the squared difference between the full spectrum and the individual slab model. The best fit (T, N, R_eff) is that which produces the smallest squared difference across all of the slab models.

As aforementioned, it is now common to fit real spectra with slab models in multiple wavelength bands (Gasman et al. 2023; Schwarz et al. 2024), or with a 2- or 3-component slab model (Temmink et al. 2024a), in order to extract information about multiple components. However, as we are anyway only modelling the contribution of warm H₂O and CO₂, we only fit a single slab, to our entire modelled wavelength range. In Appendix C.2, we test instead using the MIRI-MRS Channel 4 wavelength range (17.70–27.90 µm) for H₂O but did not find significantly different results; this illustrates how our spectra lack the full complexity of those observed (e.g. Temmink et al. 2024a).

Fig. 3 Evolution of the pebble component of the disc. Right: total dust mass (including that which has crossed the inner disc boundary and accreted onto the star) inside the initial locations of the H2O (solid) and CO2 (dashed) snow lines (right). Left: Pebble flux across the H2O snow line. Included for reference are the values from Table 6 of Romero-Mirza et al. (2024a) [RM24], which are estimated from the total number of H2O molecules <400 K, split into compact discs (small, filled, markers) and extended discs (larger, ringed, markers).

Fig. 4 Radial distribution of the H2O (blue), CO2 (bright green), and CH4 (olive green) after 1 kyr of evolution (left) and 0.1 Myr of evolution in a smooth disc in Scenario 1. Solid lines indicate molecules in the solid (ice) phase, and dashed lines molecules in the gas (vapour) phase.

4 Results

4.1 Evolution of smooth discs

We begin by briefly showing the evolution of smooth discs. To contextualise this discussion, Fig. 3 shows the evolution of the cumulative dust mass that has reached inside of H₂O and CO₂ snow lines on the right-hand panel and its rate of change, i.e. the pebble flux across the H₂O snow line, on the left-hand panel. As the CO₂ snow line is further out in the disc, there are initially more solids inside it. The pebble flux peaks at around 0.1 Myr in Scenarios 1 and 3, but somewhat later in Scenario 2. As we preserve the degree of coupling (St/α) between Scenarios 1 and 2, the Stokes number is smaller in the latter case, and the pebbles are delivered more slowly, over a longer period.

4.1.1 Spatial and temporal evolution of molecules

Figure 4 shows the radial profile in the ice and gas phase of three major molecules – H₂O, CO₂ and CH₄ – at two different times with Scenario 1 for the gas and dust transport. We see that water, being a polar molecule with a high binding energy, is only in the gas phase within around 0.47 au of the star. CO₂ is next to freeze-out, at ~2.2 au. Finally, CH₄ is the most volatile and is only frozen out outside of 16 au. After only 1000 years (left panel), there has been little redistribution of the molecules. After 0.1 Myr, when the pebble flux peaks, drift is possible for grains out to around 80 au i.e. across essentially the whole disc, resulting in a reduction in the ice abundances represented by the solid lines. For example, an amount of CH₄ ice has been moved from 16–80 au to the snow line, but it has not had enough time to spread further inwards and is thus concentrated in a small, narrow peak there. In contrast, CO₂ has been moved from a slightly wider area (2–80 au) leading to a larger peak, and the local viscous timescale is shorter at its snow line than for CH₄, and so it has spread further in. H₂O has been depleted from the largest area and so shows the greatest enhancement, and the peak has already spread so much as to produce a uniform abundance inside the snow line.

Figure 5 translates such profiles into the total vapour mass of different molecules as a function of time. Each sees an enhancement as it is delivered to its snow line by radial drift, followed by a decline once drift delivery slows and the deposited vapour is advected onto the central star. The molecules undergo this enhancement in order of their volatility. This happens because the less volatile the molecule, the smaller its initial vapour reservoir but the larger its initial ice reservoir and thus it can undergo a larger and more rapid enhancement. Moreover, pebbles start to drift at R ≳ 0.47 au and deliver H₂O to its snow line after a fairly short growth timescale, whereas dust at larger radii takes longer to grow to drifting sizes (Birnstiel et al. 2012, Eq. (13)), and so delivery to more distant snow lines starts after a short delay.

Finally, the decline from the peak enhancement is driven by the accretion of each molecule onto the star: the timescale for this to happens for molecule X can be approximated as $$t_{\text{accrete }}:=M_{\text{vap},\text{X}}/{\dot{M}}_{\text{acc},\text{X}}$$(16) $$=\frac{2\pi {\text{Σ}}_{\text{snow }}{R}_{\text{snow }}^2{A}_{\text{X}}}{3\pi {\text{Σ}}_{\text{snow }}{v}_{\text{snow }}{A}_{\text{X}}}$$(17) $$=\frac{2}{3}\frac{R_{\text{snow }}^2}{v_{\text{snow }}}$$(18) $$=0.25\text{\hspace{0.17em}Myr\hspace{0.17em}}{R}_{\text{snow }}{\left(\alpha /{10}^{-3}\right)}^{-1},$$(19)

for our choice of stellar mass and disc temperature profile, where A_X is the (assumed constant due to diffusion) abundance of molecule X in the inner disc, R_snow is the radius of the snow line in AU and Σ_snow and v_snow are the gas surface density and viscosity respectively, both evaluated at the snow line. This is of the order of the viscous timescale at the molecule’s snow line. Since the viscous timescale for advection increases linearly with radius, then the vapour of a less volatile element with a closer snow line is drained onto the star more quickly, thus also ending the period of enhancement sooner.

We thus see a sequence where the disc starts off H₂O- dominated, evolves through a CO₂-dominated phase, and finally ends dominated by the most volatile molecules such as CH₄ and CO. According to this sequence, the C/O ratio will also broadly increase over the disc’s evolution. Moreover, the timescale to go through this sequence is shorter for lower-mass stars which host colder discs with closer-in snow lines (Mah et al. 2023). Equation (19) also predicts that the time by which H₂O is removed by advection and the disc switches to the CO₂-rich phase will be later – and the subsequent lifetime of the CO₂-rich phase will be longer – at lower viscosities.

Fig. 5 Evolution of the vapour mass of H2O (blue), CO2 (bright green), and CH4 (olive green) normalised to their initial values for a smooth disc in Scenario 1. The vertical dotted lines mark the times that each molecule reaches its peak vapour mass and the rightward arrows indicate the timescale for the peak to then decrease (Eq. (19)).

4.1.2 Modelled spectrum and fits

In Fig. 6 we show an example of the method we outlined in Sect. 3, applied at an age of 10⁴ yr to the same Scenario 1 disc. Appendix C.1 discusses examples at later times of 0.1 Myr & 1 Myr, shown in Figs. C.1 & C.2 respectively.

In the top right panel, we show that the visible column densities of H₂O and CO₂ (dashed lines) as a function of radius are much lower than the total column densities (solid lines), in this case by around four orders of magnitude, as the inner disc is very optically thick and most of the molecular reservoir is therefore hidden. Moreover, while the total H₂O column density has begun to show a peak inside of the H₂O snow line (vertical blue line), reflecting the action of drift delivery as seen in the abundances in Fig. 4, this feature is not apparent in the visible H₂O column density, as the dust that delivered it contributes to increasing the local optical depth. This occurs in this scenario because the dust fragments, encounters a traffic jam¹, and thus moves at the same rate as the liberated vapour, resulting in the cancellation of the effects of enhancement and increased optical depth. The increased dust obscuration is also seen in the visible CO₂ column density, which displays a dip inside of the water snow line. Conversely, the increase in CO₂ inside its snow line is reflected in the visible column density as there is no correlated local pile up of dust.

The central panel shows channels 3B and 3C of the MIRI- MRS spectra we construct for each molecule from the temperature profile (top left) and visible column densities (top right). This particular case presents here as a CO₂-dominated spectrum due to its lack of a modelled hot water component and its fairly large CO₂ emitting area and relatively cool water emitting temperature compared to observations. The bottom panels show the χ² maps for the slab fits, with brighter colours indicating a better fit. We then overplot the parameters of the single best fit amongst these models – T, N, and R_eff – on the temperature profile and the column density profiles. For H₂O, the fit seems to trace a radius a little inside of the snow line, suggesting that it is a reasonably good measure of the area over which H₂O is in the gas phase. For CO₂ however, the radius is very similar, suggesting that CO₂ near the snow line is a bit too cold to contribute much to the spectrum. However, we do see the radius increase somewhat later in the evolution as the pile up of dust due to the traffic jam increases and the contrast between the visible CO₂ column density inside and outside the water snow line grows such that the cold CO₂ contributes relatively more to the spectrum (see e.g. Fig. C.1). For H₂O, we observe that the retrieved column density agrees well with the input visible one. For CO₂, the visible column densities are more radially variable, and the column density ends up as some sort of average. Similarly, in both cases, the retrieved temperature is somewhat higher than T(R_eff) as it is biased by the warmer gas in the inner disc. Consequently for both molecules the temperatures consistently lie in the 200–300 K range consistent with an area-weighted average $\langle T\rangle =\frac{4}{3}\frac{R_{\text{eff}}^{3/2}-R_{\text{in}}^{3/2}}{R_{\text{eff}}^2-R_{\text{in}}^2}{T}_{\text{in}}\approx 235$ K given R_in = 0.1 au, T_in = T(R_in) = 347K, and R_eff ≈ 0.32au.

4.1.3 Sensitivity of column density to underlying abundance

Figure 7 shows how the observed column densities of H₂O and CO₂ depend on their vapour masses over the course of the disc evolution, now for smooth discs in all three scenarios.

In Scenario 1, the measured water column density is almost entirely insensitive to the true vapour mass and underlying column density. As discussed above (in the context of the lack of a clear delivery signature in the visible water column density at 10⁴ yr), this occurs because both the water vapour released by the dust grains that cross the snow line – and the dust grains themselves – move with the gas inside of the snow line. Therefore, while radial drift enhances the mass of water in the inner disc, it also enhances the amount of dust and makes the disc more optically thick. The water vapour is thus hidden by the very same dust that delivers it. This constant value of the visible water column density can be explained quite simply. When the disc is highly optically thick, given that the opacity is dominated by the well-mixed small dust grains, then (see Sect. 3): $$N_{{\text{H}}_2\text{O},\text{vis}}\approx \frac{N_{{\text{H}}_2\text{O},\text{tot}}}{\tau }$$(20) $$\approx \frac{{\text{Σ}}_{{\text{H}}_2\text{O},\text{ tot }}/\left(m_{{\text{H}}_2\text{O}}\right)}{f_{\text{small }}{\text{Σ}}_{\text{dust }}{\kappa }_{\text{small }}}$$(21) $$\approx 7.2\times {10}^{19}\frac{{\text{Σ}}_{{\text{H}}_2\text{O}}/{\text{Σ}}_{\text{dust }}}{0.54}{\left(\frac{{\kappa }_{\text{small }}}{1000 {\text{cm}}^2 {\text{g}}^{-1}}\right)}^{-1}.$$(22)

Since the small grains do not evolve in size, then their opacity is fixed. Moreover, since there is a traffic jam, then dust and water vapour are transported at the same rate and hence ${\text{Σ}}_{{\text{H}}_2\text{O}}/{\text{Σ}}_{\text{dust }}$ is the same inside the snow line as when the water was in the ice phase outside the snow line. Hence, this ratio is mainly set by the initial water ice abundance on the dust grains, which we assumed to be a constant throughout the disc. As discussed in Sect. 2.3.2, observations of ices in protostellar envelopes and molecular clouds typically imply a value 2–3 times lower, while models predict one 2–3 times higher if all oxygen goes into H₂O ice. Therefore, for example, if the discs uniformly inherit the typical protostellar values, then the observed column would be lower by this factor. Nevertheless, this value can serve as a benchmark. If the column densities are much higher, it implies either a loss of inner disc dust with respect to water, or an increase in the water due to processes unrelated to dust delivery.

One way to lose inner disc dust with respect to water is if it does not couple to the gas but continues to undergo radial drift. Figure 7 shows that in Scenarios 2 and 3, where there is no change in the dust properties across the snow line and hence no traffic jam (so ${\text{Σ}}_{{\text{H}}_2\text{O}}/{\text{Σ}}_{\text{dust }}$ is no longer constant) the visible H₂O column density increases by over an order of magnitude. There are times where column density and vapour mass increase or decrease together, but during periods where the loss of accreted vapour is compensated by a loss of drifting dust (thus exposing more H₂O vapour) the result can be further regions of insensitivity to water vapour mass (Scenario 3) or even regions where vapour mass and column density anti-correlate (Scenario 2).

While Bosman et al. (2017) argued that the CO₂ emission is degenerate to changes in the global dust-to-gas ratio or CO₂ abundance, here CO₂ shows a cleaner dependence of its retrieved column density on its underlying vapour mass in all cases, though again there is no simple one-to-one relationship. This is because of the different dust dynamics – and hence evolution of the dust-to-gas ratio – at radii inside and outside the H₂O snow line. Initially, although the vapour mass starts to rise, mainly just inside the CO₂ snow line, the delivery of dust obscures the warmer and brighter CO₂ inside the H₂O snow line and the visible column there reduces a little, as seen in the above example of the spectrum construction at 0.01 Myr. As the obscured material is hotter than the delivered material, it contributes more to the line fluxes and thus the retrieved column densities are decreased overall. This effect is strongest in Scenario 1 where the traffic jam causes the dust to pile up more. However, once supply of drifting dust drops, the dust remaining between the H₂O and CO₂ snow lines drifts inwards, making this region less optically thick until it dominates the emission. Hence, the retrieved column density rises by at least two orders of magnitude, doing so more rapidly the greater the difference between the dust drift timescale and the gas accretion timescale. Finally, as the CO₂ is eventually removed by advection onto the star, the retrieved column density also drops.

We conclude that CO₂ is a somewhat more promising tracer of the enhancement by drift on timescales of ≳1 Myr since it is not so susceptible to observability effects caused by the potential coupling of the dust to the gas inside the H₂O snow line. However, the parameters of the dust and gas evolution still introduce a lot of uncertainty in the relationship between an observed column density and the underlying vapour mass.

Fig. 6 Example of the construction of a synthetic spectrum and the retrieval of parameters using 0D LTE slab model fits. The top panels show the temperature profile (left) and the column density profiles after 104 yr (right) for small dust (black) and total (solid) and visible (dashed) H2O (blue) and CO2 (green). The light, vertical, dotted lines indicate the snow lines. The middle panel shows the synthetic spectra of each molecule (assuming a distance of 140 pc) while the bottom panels show the goodness of fit for a grid of slab models with varying temperature and column density. The white contours show the radius (in AU) of a circle whose area maximises the goodness of fit for each combination of temperature and column density. The overall best fit parameters are shown on the top panels as the crosses, with the vertical error bars representing the spacing of the grid of slab models.

Fig. 7 Relationships between the vapour mass and the column density measured using slab fits for H2O (left) and CO2 (right) for a smooth disc model in Scenarios 1 (purple), 2 (gold), and 3 (red). The curves shown are traversed anticlockwise.

4.2 Evolution of gapped discs

4.2.1 Scenario 1: Discs with traffic jams

We now present the full results for our first parameter set, which leads to the traffic jam inside the snow line. Figure 8 shows the evolution of the vapour masses and retrieved column densities for models with different gap radii. As discussed, in the smooth disc, the H₂O vapour mass rises more quickly than that of the CO₂. This is reflected in an initial drop in the CO₂/H₂O vapour mass ratio, followed by a subsequent increase as the H₂O vapour mass starts to decline. At the latest times, the models tend back towards the initial ratio again as then the main source of both molecules is sublimation from remnant small grains that are fully entrained with the gas and so do not deliver H₂O to the inner disc any faster than the gas delivers CO₂.

Including a gap in the disc cuts off the pebble flux once the growth front has reached the gap location, thus reducing the delivery of molecules to the inner disc. Except for the largest R_gap cases, this happens before the peak enhancements of H₂O and CO₂ expected in a smooth disc, which lowers the maximum enhancement and shifts the peak earlier. Generally speaking, the closer in the trap, the sooner it interrupts the dust flux, and the stronger these effects are, as was shown for H₂O by Kalyaan et al. (2021, 2023). Moreover, as H₂O is depleted more quickly by advection than CO₂, it reacts faster to the loss of supply. As a result, traps act to increase the CO₂/H₂O vapour mass ratio, which now peaks in the range 30–300, compared to at around 10 for a smooth disc. This is consistent with Mah et al. (2024) where sufficiently deep gaps raised the inner disc C/O.

At late times, the models with traps produce vapour masses that flatten off. While smooth discs are constantly losing both gas and dust to the star, discs with traps by definition are keeping hold of a much larger reservoir of dust, which slowly leaks through the traps and thus sustains a minimum supply of vapour mass to the inner disc. We explore this more in Sect. 4.3. Generally the dust that leaks grows inefficiently and remains small and well coupled. Hence in this phase, there is no difference between the transport of vapour and ices, so CO₂ and H₂O are delivered at the same net rate, fixed by their ice abundances, and thus their ratio tends back towards its initial value.

One special case is where the trap is sited at 1 au, which is outside the H₂O snow line but inside that of the CO₂. In this case, it has little effect on the CO₂ vapour evolution but a very significant one on the H₂O vapour evolution. This results in a much earlier rise in the CO₂/H₂O vapour mass ratio, which reaches ~1000 on timescales ~0.3 Myr, higher than achieved in any of the other models by a factor of a few. The 1 au case also shows, uniquely, a resurgence of the H₂O vapour mass due to leaking dust. This happens largely because the dust mass is still accumulating in the trap on million year (Myr) timescales. This long-term supply of H₂O leaking from the trap is similar to that found for traps of intermediate depth by Mah et al. (2024, their so-called ‘traffic jam scenario’) which balance the ability of deep traps to retain some dust with the ability of shallow traps to then leak that dust. However, because all the CO₂ enters the gas phase before the pebbles are trapped, the dust in the trap is CO₂-poor. Thus, the slow leaking of dust from the trap is not able to replenish the CO₂ vapour, resulting in very low CO₂/H₂O ratios on timescales ≳3 Myr.

The column densities tell a similar story. Despite the reduced dust flux in the presence of traps, the inner disc never becomes sufficiently optically thin to deviate from Eq. (22), and thus in all cases the H₂O column density is fairly constant in time, varying by no more than factor ~3, considerably less than the variation in the underlying H₂O vapour mass. On the other hand, despite the fact that traps reduce the delivery of CO₂ to the inner disc, by blocking dust from reaching the inner disc they also allow the visible column density to undergo its rapid rise much sooner and thus tend to promote higher retrieved CO₂ column densities. Moreover, as the dust obscuration affects both the H₂O and CO₂ in a correlated way, when the CO₂/H₂O column density ratio is considered, the evolution (for any gap location) closely mirrors that of the vapour mass ratio.

With the caveat that our spectra lack contribution from the hot surface layers (see discussion in Sect. 3), we collect a sample of discs with so far published or submitted analysis of their MIRI-MRS spectra (Table 5) where slab model fits for both H₂O and CO₂ are available. We limit ourselves to discs around stars >0.3 M_⊙ because lower mass stars are known to host discs with a more hydrocarbon-rich, H₂O-poor, chemistry, which may be connected to the fact that such discs are cooler and their snow lines lie closer to the star, and hence the timescales on which their chemistry evolves due to transport are much shorter (Mah et al. 2023) than for our fiducial 1 M_sun case. Strictly speaking, given that we do not model the hot surface layer of the disc or the extended snow surfaces, we should attempt only to compare to the warm (250–550 K) H₂O and CO₂ components. However given the already small sample this leaves us, and given that the works all differ in the wavelength ranges and number of components they fit for the water (and in any case only ever fit a single CO₂ component that may or may not be within the temperature ranges we find) we simply plot an unweighted average over all of the H₂O fits in each work, and indicate their range with the error bar. By so doing, we make the assumption that the chemistry of all the components is correlated, for example that H₂O enhancement deep in the disc would also enhance H₂O in the surface layers (rather than being completely reset by reactions), and so that all components contain information about the bulk chemical composition even if they do not trace the bulk material directly.

For both H₂O and CO₂, we see that our models generally predict column densities at the upper limit of what is observed. As discussed, in the case of H₂O, the value depends simply on the initial (H₂O) ice-to-dust ratio and the opacity; the CO₂ column density will also be affected by equivalent quantities. Lowering the initial ice abundances, as might be more consistent with protostellar estimates, and/or raising the dust opacity, could potentially account for some of the difference between the models and the bulk of the observations. Alternatively, as the dominant line features in the slabs are often very optically thick (as noted for DR Tau by Temmink et al. 2024a), the observations may instead only trace the column down to where these features – rather than the continuum – become optically thick. For CO₂, better agreement is seen between the models and the more CO₂-rich discs, particularly when the ¹³CO₂ isotopologue estimates are used. This may reflect the fact that ¹³CO₂ components are colder (Vlasblom et al. 2025) and may potentially therefore originate deeper within the disc (Bosman et al. 2017, 2018, 2022b), consistent with their expected lower (line) optical depth. They may therefore be more reflective of the bulk material that we try to model rather than the warm surface chemistry. Indeed the observed temperatures of these components are closer to those that we model. Finally, the predicted column density ratios are much closer to those observed, suggesting that the retrieved column densities of the two molecules may be less than our predictions due to similar effects.

Focussing therefore on ratios, firstly for two discs with ¹³CO₂ detections – GW Lup (Grant et al. 2023) and CX Tau (Vlasblom et al. 2025) – the former is most consistent with a CO₂/H₂O column density ratio predicted by a model with a gap at large radii – consistent with the presence of a single known gap at 74 au (Huang et al. 2018) – while the latter is more consistent with a gap closer than 7 au, or a completely smooth disc – consistent with its status as a compact, ‘drift-dominated’, disc (Facchini et al. 2019) with no resolved substructures down to 5 au. The water-dominated DR Tau (Temmink et al. 2024b,a) is a compact disc, with no substructures that are clearly resolved in the image plane. Our models tentatively suggest that it is most consistent with a gap very close to the star, between the H₂O and CO₂ snow lines, that is, leaking considerable amounts of H₂O on Myr timescales. The other discs are all known to be substructured, though AS 209 has a much more complicated system of multiple rings (Huang et al. 2018) which may be why it agrees the least, and SY Cha is a transition disc with a large millimetre dust cavity (Schwarz et al. 2024). Overall, this may suggest that CO₂-dominated spectra reflect discs where transport controls the chemistry while other spectra represent scenarios where chemistry is determined by other processes such as chemical processing in the warm molecular layer. Alternatively, the lack of clearly detected hot bands in most of these discs (due to strong blended H₂O emission) may bias the observations away from high column densities. Grant et al. (2023); Gasman et al. (2023) show that the Q-branch alone can be fit with an optically thin, low column density solution, but that including the hot bands leads to column densities more than an order of magnitude higher in which the Q-branch is optically thick and insensitive to column; this would explain why Sz 98, where those bands are tentatively detected, shows a column density more similar to the discs with CO₂-dominated spectra and ¹³CO₂ detections.

We find that, especially in the case of water, the effective emitting radius and temperature derived from the slab fits are mostly determined by the location of the snow line and the average disc temperature interior and show almost no evolution (e.g. Fig. C.4). Therefore, we do not include plots of T or R_eff for our full grid of models.

Fig. 8 Evolution of the vapour mass (left) and column density (right) for CO2 (top row), H2O (middle row), and their ratio (bottom row) for discs with traffic jams inside the water snow line (Scenario 1). The grey line represents a smooth disc with no traps, while the blues lines, in order of increasingly dark shade, represent traps with increasingly close in traps. Overplotted at their approximate age are the results of slab model fits for a few discs from the literature: CX Tau (Vlasblom et al. 2025), GW Lup (Grant et al. 2023), SY Cha (Schwarz et al. 2024), Sz 98 (Gasman et al. 2023), DR Tau (Temmink et al. 2024b,a), AS 209 (Romero-Mirza et al. 2024b) and DoAr 33 (Colmenares et al. 2024), ordered by increasing stellar mass (more yellow colours). The open squares indicate an alternative estimate obtained from the 13CO2 isotopologue. The grey bar in the top left shows the ± log10(2.5) uncertainty depending on the basis used for the initial ice abundances (Sect. 2.3.2).

4.2.2 Scenario 2: Rapid drift and slow gas accretion

In Scenario 2, there is no traffic jam inside of the snow line and therefore no inevitable hiding of the water. Thus, once the dust flux drops, the dust inside the snow line also drifts onto the star, and the delivered water is rapidly revealed. At the same time, the gas is no longer being so effectively advected onto the star, and so the vapour masses build up much more, peaking later. The result of the loss of opacity, revealing deeper layers of the disc, is an almost vertical climb in the column densities within 1 Myr to values $N_{{\text{H}}_2\text{O}}\gtrsim {10}^{24} {\text{cm}}^{-2}$ and $N_{{\text{CO}}_2}\gtrsim {10}^{22} {\text{cm}}^{-2}$, far in excess of the observed ranges for discs at ~1 Myr. The slower gas evolution means that the general pattern undergone by the CO₂/H₂O ratio (in both vapour mass and column density) of an initial fall and subsequent increase (once H₂O is accreted) happens much later (but otherwise looks similar). Therefore, by 1–3 Myr, all models are still recovering from the initial drop in the H₂O/CO₂ ratio as there has been less removal of the water by accretion onto the star. Thus, the ratio of the retrieved columns is brought much closer to the low values observed for CO₂-poor discs around AS 209 and SY Cha. However, the high absolute columns, and the inability to replicate the more CO₂ -rich sources, suggest this scenario is unlikely to be the correct explanation without invoking very efficient chemical reprocessing and low resulting molecular abundances.

4.2.3 Scenario 3: Rapid drift and faster gas accretion

Much like in Scenario 2, in Scenario 3 there is no traffic jam inside the water snow line and therefore the column densities increase sharply after the pebble flux is interrupted by the trap. However the assumptions of different α values for turbulence and viscosity means that the vapour mass enhancements do not build up quite so high or last quite so long, more like Scenario 1. The result is somewhat more reasonable column densities of 10²⁰ −10²¹ cm⁻² on Myr timescales, though these are still a bit higher than observed values, similarly to the first scenario. Overall, this scenario and the first perform similarly well at reproducing the CO₂/H₂O ratio in ~1 Myr discs, and we can draw the same conclusions about the preferred locations of any gap amongst our models. However, in this scenario there is more sensitivity to variations in the H₂O abundance than in Scenario 1.

4.3 Trends with gap location

As could be seen above, the location of any gap in the disc affects the evolution of the vapour masses and column densities of different molecules through its effects on the time that the pebble flux is interrupted, the amount of pebbles that are thus isolated from the inner disc, and the amount of this isolated material that leaks through the trap. We also suggested that between two CO₂-dominated sources observed with JWST, GW Lup seems to prefer a more distant gap location and CX Tau a closer gap location (or none at all). In this section we therefore explore in more detail the relationship between gap location and the vapour masses or column densities, focusing on Scenario 1 which produced the most reasonable column densities (Scenarios 2 and 3 can be found in Appendix B). Figure 9 shows that for the vapour masses of each molecule, there are clear trends in the vapour masses with radius. These can be divided broadly into positive relationships with R_gap at early times in the disc evolution and negative relationships at late times.

Fig. 9 Trends of CO2 (top) and H2O (middle) vapour masses (left) and column densities (right) and their ratio (bottom) with Rgap at different stages of evolution for discs with traffic jams inside the water snow line (Scenario 1). The initial condition is in grey, while increasingly dark shades of green correspond to increasingly late times. The model with Rgap = 1 au is plotted separately at the dots as its gap is located very close to the H2O snow line and thus inside the CO2 snow line. A smooth disc model is included for comparison as the triangles. Equations (23) (dashed-dotted) and (11) (dotted) are plotted as guides to the eye the trends for approximate trends at early and late times.

4.3.1 Trends in vapour mass at early times

The positive relationship between the vapour mass and gap location at early times ≲0.5 Myr, which is approximately linear, can be easily understood as tracing the amount of mass that initially lies inside of the trap and is thus transported to the inner disc unimpeded: $$\begin{array}{l}{M}_{\text{d},\text{in }}=M_{\text{d},\text{tot }}\left(1-\exp \left(-\frac{R_{\text{gap}}}{R_{\text{C}}}\right)\right)\\ \approx M_{\text{d},\text{tot }}\frac{R_{\text{gap }}}{R_{\text{C}}}\left(R_{\text{gap }}\ll R_{\text{C}}\right).\end{array}$$(23)

The dependence predicted by Eq. (23) is plotted in Fig. 9 (dashed-dotted line) and is indeed roughly consistent with the models, especially for CO₂. However, the H₂O vapour delivered to the inner disc accretes on timescales ∼0.1 Myr, several of which have passed by the first snapshot at 0.5 Myr. Hence in these cases, the relationship is already evolving towards its late time behaviour, making it deviate more significantly from the approximately linear trend, though it nevertheless retains a positive trend up until a time of ∼1 Myr. Thus, the closer to the star the gap is opened, the less dust there is inside of it to drift, and thus we get less enhancement of the vapour masses, as noted for mm-sized pebbles by Kalyaan et al. (2021).

4.3.2 Trends in vapour mass at late times

Once the molecular reservoir that was originally present inside the gap has all been accreted onto the star, both the CO₂ and H₂O vapour masses transition to a negative dependence on R_gap . At these late times, as mentioned briefly in terms of the flat vapour mass evolution, the vapour in the inner disc is supplied by whatever small dust manages to leak through the trap. Traps nearer to the star are leakier (Kalyaan et al. 2023) and so more dust is escaping.

Although small dust dynamically mostly follows the gas, it can also be produced by the fragmentation of the large grains that are tightly concentrated in the trap. Although we do not model this directly, the fixed f_small means that we will also have a population of small dust that follows the concentration of the large dust grains. This results in diffusion being a strong driver of small dust leaking through the trap. Using Eq. (11), we can estimate a diffusive flux of $$F_{\text{leak }}\approx 2\pi R_{\text{gap}}\frac{v}{Sc}\text{Σ}\frac{{ϵ}_{\text{small}}}{w_{\text{d}}}$$(24) $$={\frac{2{\dot{M}}_{\text{acc }}{f}_{\text{small }}}{3Sc}\frac{R_{\text{gap }}}{w_{\text{d}}}\frac{{\text{Σ}}_{\text{d}}}{{\text{Σ}}_{\text{g}}}|}_{\text{gap }},$$(25)

where є_small = f_smallє is the mass fraction in small dust for є the total dust mass fraction, and where we use the steady state relationship Ṁ_acc = 3πνΣ_G. w_d is the radial width of the concentrated dust in the trap (set by the radial width of the fragmenting large grains), which can be related to the width of the pressure bump in gas by $w_{\text{d}}^2=\frac{\alpha }{St}{w}_{\text{g}}^2$ (e.g. Rosotti et al. 2020); we shall further assume w_g ≈ w_gap . For the dust-to-gas ratio, we consider M_d,out = M_d,tot − M_d,in (Eq. (23)) to be concentrated in a ring of area 2πR_gapw_d, which gives ${\frac{{\text{Σ}}_{\text{d}}}{{\text{Σ}}_{\text{g}}}|}_{\text{gap}}\approx {ϵ}_0\frac{R_C}{w_{\text{d}}}\exp \left(-\frac{R_{\text{gap}}}{R_{\text{C}}}\right)$ where ϵ₀ ≈ 0.01 is the initial dust mass fraction.

The growth of large grains in the trap will be limited by turbulent fragmentation, in which case Eq. (6) predicts that $\frac{\alpha }{St}\propto c_{\text{S}}^2\propto R_{\text{gap}}^{1/2}$ ,. Combining all of the above gives $$F_{\text{leak }}\propto \frac{R_{\text{gap }}}{w_{\text{d}}^2}\exp \left(-\frac{R_{\text{gap }}}{R_{\text{C}}}\right)$$(26) $$\propto \frac{R_{\text{gap }}^{3/2}}{w_{\text{gap }}^2}\exp \left(-\frac{R_{\text{gap }}}{R_{\text{C}}}\right),$$(27)

where we assume f_small to be constant as per the two-population model (Birnstiel et al. 2012)².

Therefore, since our models assume $w_{\text{gap }}\propto H_{\text{gap }}\propto R_{\text{gap }}^{5/4}$, then this gives F_leak ∝ R⁻¹ for R_gap ≪ R_C when turbulent fragmentation limits dust growth throughout the trap. The limit from Eq. (27) is plotted as the dotted line in Fig. 9 While the agreement is imperfect as (a) we made several approximations in the calculation of the dust flux and (b) the two population model does not completely accurately capture the nuances of diffusive flux, there is sufficient agreement with the general negative trend to be confident that diffusion of dust from the trap is the explanation for the trend.

In conclusion, the negative trend with R_gap results principally from the trapped dust being more tightly concentrated in traps closer to the star, which induces a stronger diffusive gradient.

4.3.3 Observable trends

As before, the H₂O vapour evolution leads that of the CO₂. This can be seen most clearly on timescales of a couple of Myr, where H₂O has adopted a mostly negative relationship with R_gap (save for the cases with the most distant gaps where the most volatiles are delivered and consequently it takes the longest to deplete them to the same level) while the CO₂ still has a positive trend. This is important because it means that the CO₂/H₂O ratio is generally an increasing function of R_gap (though declines again slightly for R_gap ≳ 60 au as the gap location becomes close to the disc outer edge at R_C = 70 au) on timescales of 1–3 Myr. However, at late times, the ratio is generally independent of trap radius (except for if the trap lies between the two snow lines) and is simply determined by the input ice abundance ratio because the supply of molecules is controlled by leaking from the traps and hence both H₂O and CO₂ leak at the same rate.

If measurable, the existence of distinctly opposite trends at the early and late stages of disc evolution would be a powerful one. It would potentially allow us to measure the age of a population of structured discs and their traps relative to the timescales of the dust and gas evolution that control the trends. Therefore, we finally assess whether such trends carry over to the column densities of the molecules retrieved from slab fits.

As expected (given the previously discussed insensitivity of the column density to the vapour mass in Scenario 1), the trends with respect to H₂O are lost. Conversely, we generally recover the positive relationship expected between the CO₂ column density and gap radius at early times, as there is no dust pile up between the CO₂ and H₂O snow lines, giving us sensitivity to the effect of close in traps reducing the enhancement of CO₂ by dust delivery. However, at late times, we see the same column densities emerging for all traps, despite the fact that the closer in traps are leakier. This is because the leakier traps are supplying more small dust – all of which remains coupled to the gas – which contributes to increasing the continuum optical depth and obscuring the delivered CO₂ (much like applies inside the H₂O snow line at all times due to the traffic jam). Finally, because the H₂O column density never shows any trend, these trends for the CO₂ column density carry over to the ratio of the CO₂ and H₂O column densities. Given that dust optical depth ultimately affects both molecules similarly, the trends in column density ratio reflect the underlying CO₂/H₂O vapour mass ratio relatively well.

Fig. 10 Dependence of the CO2/H2O column density ratio on gap location at different snapshots in the disc evolution. The panels show models in which the traps are fully developed at successively later times.

4.3.4 Dependence on trap formation time

Having identified the CO₂/H₂O column density ratio as a reasonably promising tracer of the role of drift and traps preventing drift, we now consider the effect of traps forming later (up to 1 Myr into the disc’s evolution) on its predicted trends with R_gap.

The time taken for dust to start drifting is generally an increasing function of radius, as the dust is assumed to start small and need some time to grow, with the growth timescale which delays the onset of drift scaling as 1 /Ω_K ∝ R^3/2 (Birnstiel et al. 2012). As long as a trap forms before the dust at that radius (and exterior) starts to drift, then it will have the same effect as a trap formed at t = 0. However, if the dust has already started to drift past that radius once the trap forms, then the amount of dust it can prevent reaching the inner disc is reduced. Combining these considerations, we can see that delaying the onset of trap formation will have the strongest effect on traps closer to the star where dust starts to drift sooner.

The trend at early times is controlled by the amount of dust that is not blocked by the trap, which is reduced most by closein traps. Later trap formation therefore allows more dust into the inner disc before close-in traps block it, leading to more enhancement of the inner disc; the outer traps are not so affected because most dust lies interior to them anyway. The result, as illustrated by Fig. 10, is to make all gaps act more like the gaps further out in terms of their impact on CO₂/H₂O. This makes the positive trends with R_gap shallower: if the traps still form within ≲0.1 Myr (the rough growth+drift timescale for the outer disc) a majority of the dust will still be outside any trap once it opens, and a positive, albeit less strong, relationship with R_gap on ∼2 Myr timescales is still maintained. On the other hand, such trends are either largely indiscernible or even reverse if the traps form after ≳0.1 Myr. Likewise, Mah et al. (2024) found 0.1 Myr to be a critical timescale for trap formation to prevent significant H₂O enhancement, and showed that otherwise the enhancement would last the whole disc lifetime for α ≲ 10⁻⁴.

There are known examples of rings suggestive of dust traps in Class 0/I discs at young ages (ALMA Partnership et al. 2015; Segura-Cox et al. 2020), but such cases are rare, with surveys suggesting only 10% incidence for substructures ≳10 au (Ohashi et al. 2023; Hsieh et al. 2024). If any connection can be demonstrated between dust substructure and inner disc chemistry, this would imply a significant population of unresolved or unobservable substructures in these discs, potentially due to the dust emission being both optically and geometrically thick (Ohashi et al. 2023; Sharma et al. 2023; Guerra-Alvarado et al. 2024). This would be consistent with the population synthesis modelling of Delussu et al. (2024) which requires substructures capable of trapping dust to be nearly ubiquitous and to form by 0.4Myr, in order to match the millimetre spectral indices, fluxes, and sizes of discs. However, it is hard for planets to grow large enough to open gaps at such times, especially given that younger discs are warmer; even with the generous assumption that planetary embryos can form in the Class 0 phase, less than 10% of discs are likely massive enough to produce gap-opening planets in under 0.1 Myr (Nazari et al. 2024). Indeed, although the relative mass budgets of Class II discs compared to planetary systems (Manara et al. 2018) and embedded discs (Tychoniec et al. 2020) point towards the solids that become planets being locked up into larger bodies before the Class II stage, this may only have time to progress as far as planetesimals in the embedded phase (Cridland et al. 2022). Consequently, a strong connection between dust substructure and inner disc chemistry would likely require a non-planetary origin for the substructures such as magnetohydrodynamic zonal flows (see discussion in Bae et al. 2023), or instabilities related to magnetohydrodynamic winds (Riols & Lesur 2019).

5 Discussion

5.1 What causes CO₂-dominated spectra?

Our models show that the inner disc CO₂ /H₂O ratio in the bulk of the disc can be traced through 0D LTE slab model fits to MIR spectra with components of 200–300 K under assumptions of efficient vertical mixing and negligible chemical reprocessing. This ratio is less affected by continuum opacity effects than the individual columns and therefore is potentially a better tracer of the inner disc chemical evolution as a result of gas and grain transport than the features of the individual molecules. This is because the abundance of the two molecules evolves on different timescales, leading to a predictable pattern in the ratio of their vapour masses due to drift while the relative evolution of the molecular gas and the dust is more uncertain. The column density ratio closely reflects the pattern shown in the vapour masses for a range of dust evolution scenarios.

5.1.1 Measuring CO₂/H₂O and Limitations of 1D models

To use the CO₂/H₂O ratio as a tracer of drift, one needs to be able to retrieve it reliably from observed spectra. Our models, at the 1–3 Myr ages of most discs characterised so far, mostly predict higher column densities than observed; the predicted temperatures are also lower than the observed temperatures.

One possible reason for this would be that the observed spectra trace a relatively thin vertical layer that is dominated by the disc’s surface layers; here gas temperatures are hotter and molecules are exposed to higher UV fluxes, leading to faster chemical destruction and a reduction in molecular abundance. In other words, the disc surface has much shorter chemical timescales than estimated at the midplane in Appendix A, whereas our model ignores chemical reactions entirely. This would be in line with the conclusions of Bosman et al. (2018) that a vertically averaged destruction rate at least 1000 times higher than that at the midplane would be needed to explain the inferred CO₂ abundances.

Nevertheless, the existence of relationships between the cold H₂O emission and the disc structure (Banzatti et al. 2023; Romero-Mirza et al. 2024a) – while the hotter emission component shows no such relationship but instead correlates with the accretion rate – suggests that chemical reprocessing in the surface layers does not result in the complete destruction of a drift signal with some signature remaining in the colder emission. Even if the exact $N_{{\text{CO}}_2}/N_{{\text{H}}_2\text{O}}$ is altered, it should still correlate with C/O, which evolves following the delivery and accretion of H₂O and CO₂ (Mah et al. 2023). Put another way, any signature of delivered CO₂ should last longer than any signature of delivered H₂O. Nevertheless, 2D thermochemical modelling is needed to determine how the delivered ratio translates into the observable ratio and make the probe more quantitatively robust.

A second difference between our modelling approach and real spectra is that we always use the entirety of the molecular spectrum within our modelled wavelength range. For example the full structure consisting of several bands and branches of CO₂ emission shown in Figs. 6, C.1, and C.2 may allow us to break the degeneracy between column density and temperature for optically thick features (Bosman et al. 2017, 2022b). Conversely, in many observed spectra only the 15 µm Q-branch of CO₂ is unambiguously identifiable amongst the forest of H₂O lines.

Notably, the discs that come closest to our models are those where additional CO₂ features are seen. Grant et al. (2023) and Vlasblom et al. (2025) show that the estimated column densities are much higher when the weaker ¹³CO₂ Q-branch feature is used and imply the main isotopologue has become optically thick. This is consistent with modelling predictions that the main contribution to the ¹³CO₂ emission should come from slightly deeper in the disc (Bosman et al. 2017, 2018) and be a sensitive tracer of drift, a hypothesis that is also supported by the retrieved excitation temperatures of ¹³CO₂ being colder than those of CO₂, more similar to our modelled temperatures.

The degeneracy between the effects of column high optical depth and high temperature on the broadening of the CO₂ 15 µm Q-branch means that this feature could be fit with a hot, optically thin, solution (Salyk et al. 2011). However, Grant et al. (2023) demonstrate how including the hot bands lowers the retrieved temperature and raises the retrieved column density substantially. A similar effect was seen by Gasman et al. (2023), and Colmenares et al. (2024) also include a hot band in their fits and find an above average $N_{{\text{CO}}_2}/N_{{\text{H}}_2\text{O}}$. This is consistent with the expectations from thermochemical modelling that the hot bands trace deeper emission than the 15 µm Q-branch (Bosman et al. 2022b). Finally, for fits to protostellar spectra, the P and R-branches alone give a much lower temperature and higher column than when the Q-branch is also included (Fig. D.1; van Gelder et al. 2024).

Overall, it appears key that weaker features sensitive to higher molecular column densities and deeper depths in the disc are used to mitigate against the retrieval being dominated by hot, chemically altered gas at the disc surface and thus estimate the bulk ratio more accurately. For CO₂, this is likely mainly possible for already quite CO₂ -rich discs where these features are clear amongst the forest of other emission lines and can more easily dominate the spectrum.

5.1.2 Interpreting CO₂/H₂O

In smooth discs, both H₂O and CO₂ are delivered efficiently by dust to the inner disc, but water is accreted onto the star faster. Therefore, the CO₂/H₂O ratio tends to be enhanced by radial drift on timescales 0.5–1 Myr (Mah et al. 2023) before declining back towards a baseline. Thus, even smooth, drift-dominated discs can produce CO₂-dominated spectra when observed on Myr timescales (assuming a solar mass star and α ~ 10⁻³), as shown by the reasonable agreement between our model and the observations of CX Tau (Vlasblom et al. 2025) when using the estimate from the ¹³CO₂ isotopologue. In this picture of sequential delivery as a result of drift, more H₂O-rich spectra may represent discs that are younger or undergoing slower evolution (e.g. Scenario 2). Relatively more evolved sources should show weaker H₂O emission while the CO₂ has not yet declined, allowing clearer identification of weak CO₂ features that typify CO₂-dominated discs.

We also found that dust trapping tends to enhance CO₂/H₂O. As the H₂O is accreted onto the star faster than CO₂, the best way to achieve this on Myr timescales is to have interrupted the delivery of H₂O early enough for it to have all accreted, but late enough to still allow plenty of CO₂ to have been delivered. Consequently, traps at 10s au (which still let the early drift phase happen but do not leak much dust and ice at later times and therefore cannot replenish the accreted water) are more effective at enhancing CO₂/H₂O than traps at small radii. However, one special exception is where the trap lies between the snow lines of H₂O and CO₂ in which case CO₂ arrives unimpeded to the inner disc, while almost no dust can reach the water snow line.

Different discs also show variations in the temperatures and emitting areas which potentially point towards a difference in the 2D thermochemical structures of their inner regions. One way in which the inner disc structure can be altered such as to affect the MIR spectrum is with to the presence of a deep gas cavity, unresolved in the ALMA observations (Vlasblom et al. 2024). While small cavities can increase the emitting area due to the heated inner wall, cavities of a few AU in size reduce the area over which the disc is hot enough to desorb ices and also form molecules through gas phase reactions; this tends to enhance the CO₂ line strengths relative to those of H₂O. However, the impact on retrievable column densities was not explored, and these may be more weakly affected if the visible column density is a weak function of radius as in Fig. 6. While cavities nevertheless remain a viable option, our models suggest an alternative route to producing CO₂-dominated spectra; future high-spatial- resolution submillimetre observations would help to distinguish the two scenarios.

One mechanism that may lead to the formation of cavities of a few au is internal photoevaporation, with this occurring once the accretion rate falls below the photoevaporation rate (Clarke et al. 2001). In such a case, the inner disc is completely isolated and thus no leaking of dust into the inner disc is expected once the cavity opens. Lienert et al. (2024) explored this scenario in a disc undergoing X-ray photoevaporation (Picogna et al. 2021) and with α = 10⁻⁴ and hence slow viscous evolution (c.f. our Scenario 2). In this case, they find that photoevaporation opens a cavity before the CO₂-rich vapour reaches the inner disc (a few Myr, see Fig. B.1) and hence the CO₂-dominated phase is prevented. Moreover, given the weak accretion implied by a low α, the sublimated water is trapped in the vicinity of the snow line by the cold finger effect and is poorly accreted, leading to a resurgence of the H₂O abundance in the inner disc during the accretion of the inner disc. Thus, whether cavities indeed lead to CO₂-dominated or H₂O-dominated spectra depends sensitively on the timing of their formation and the relative roles of equilibrium (photo)chemistry and transport.

Finally, Vlasblom et al. (2024) also show that having more small dust present in the upper layers reduces the temperatures and therefore the size of the hot surface layer and the abundance of H₂O relative to CO₂ therein. This could also lead to CO₂- dominated spectra. Such dust enrichment could also be the result of drift (as suggested for IM Lup by Bosman et al. 2023), and so this effect may broadly act in the same direction as our results. Future modelling should therefore consider if the temperature structure of the inner disc changes as a result of dust drift.

5.2 Survival of dust upon ice sublimation

The outcomes of our modelling, particularly the sensitivity to the enhancement of H₂O by radial drift (as well as the late time behaviour dominated by the leaking of small dust from the traps), are dependent on how well-coupled the dust grains are in the inner disc. In the two-population model (Birnstiel et al. 2012) on which we base our dust, it is assumed that dust can grow to an equilibrium large size resulting from coagulation and – unless drift removes grains faster – collisional fragmentation. In Scenario 1, this leads to a change in the dust size across the H₂O snow line due to the assumption that dry grains are more fragile. However such an approach neglects the fact that dust grains can, realistically, have a much more complicated internal structure: for example if they coagulate while completely covered in ice, then the sublimation of this icy glue could lead to the dust grains fragmenting. This is the ‘many-seeds’ model (Schoonenberg & Ormel 2017), which may be expected to produce well-coupled micron-sized grains (Saito & Sirono 2011), much smaller than the fragmentation limit. Therefore, even in situations (e.g. Scenarios 2 and 3) where fragmentation-limited dust grains drift, the ‘many-seeds’ model would predict that unless the dust grains can quickly regrow to large sizes (which we see for the small dust leaking from the traps in Scenario 3 but not in Scenarios 1 or 2, Appendix B); the outcome could be more like Scenario 1, where the dust and vapour have the same evolution and we would lack sensitivity to the evolution of the H₂O.

In this context, Houge & Krijt (2023); Houge et al. (2024) have sought to differentiate these scenarios through modelling the collisional evolution of dust in outbursting discs, where the disc undergoes transient heating, temporarily leading to sublimation over a larger area. Although Houge & Krijt (2023) argue that regrowth to equilibrium may take only several thousand years to re-establish at a few au, which is comparatively short compared to the accretion timescale, the low continuum spectral indices of the disc of the currently outbursting V883 Ori imply the presence of cm-sized grains and thus require that dust grains be indeed resilient to breakup during sublimation (Houge et al. 2024), meaning that the many-seeds model does not apply.

5.3 Priorities for future studies

In order to properly test the scenarios put forward herein of the sequential delivery of volatiles to the inner disc following drift and the impact of dust traps on this picture, we suggest that observations spanning a greater range of disc ages are necessary. In our fiducial Scenario 1 at least, the period of 1–3 Myr where most data so far have been taken sees a rapid decrease in the CO₂/H₂O ratio – the best tracer that we find of the true evolution – resulting in models spanning a large range of values in line with a large scatter in the data, making it hard to discern any real evolution. Observations around the expected peak of this ratio before 1 Myr, as well as at ages older than 5 Myr would much better constrain the chemical evolution of the disc. Multiple accepted Cycle 2 JWST programmes (e.g. GO 2970- Pascucci et al. 2023, GO 3034- Zhang et al. 2023) should prove useful with regards to the older discs. Moreover, observations of even younger discs would also be important to test whether discs do mostly go through a water-dominated phase earlier on.

Secondly, the comparison of models to data is made challenging by the lack of consistency as the analysis methods have rapidly developed (e.g. Temmink et al. 2024a; Kaeufer et al. 2024) since the first JWST disc spectrum published (Grant et al. 2023). In particular, different works differ in how column densities are extracted from those that are available in terms of the number of slabs that are fit, or the spectral ranges or features to which they are fit. In order to argue for an underlying cause of a difference, for example that the higher CO₂/H₂O ratio in GW Lup than CX Tau reflects the action of a dust trap at large radius compared to that of a trap at a small, unresolved radius (or no trap at all), one needs to ensure that the values are derived consistently, and to understand better what part of the molecular budget of the disc the measurements probe. In order to refute or confirm the predicted relationships between gap location and CO₂/H₂O and therefore the role of drift and trapping in controlling inner disc chemistry, we conclude that a priority on the observational side should be to conduct a more systematic, coherent, sample analysis of JWST spectra at different phases of protoplanetary disc evolution using the weaker features that best represent the bulk material below the hot disc surface layers, while 2D thermochemical modelling needs to demonstrate the spatial origin of different spectral components and help justify the spectral analysis. Moreover, modelling should consider the influence of the evolving midplane chemistry on the more easily traced surface layers as a result of vertical transport and subsequent reactions, to allow the column density ratio to be used more quantitatively as a tracer.

6 Conclusions

Previous works have explored the link between dust drift and traps on the abundance of H₂O (Kalyaan et al. 2021, 2023) and the overall C/O ratio (Mah et al. 2023, 2024) in the inner regions of protoplanetary discs. We have built on these studies in order to understand how these processes also affect the abundance of CO₂ in these inner regions, and, crucially, how strongly these molecules would appear in MIR spectra depending on the disc’s age and the time and location of trap formation. Our main findings can be summarised as follows:

1. Ices sublimating from drifting dust as it crosses snow lines, followed by the gaseous advection of the vapour inwards, lead to the enhancement of the inner regions of protoplane- tary discs with molecules. This is followed by the accretion of the molecule-rich gas onto the star. This sequence proceeds faster for less volatile molecules, leading to a H₂O-rich stage followed by a CO₂-rich stage (Fig. 5).

2. How well this evolution of the molecular content is traced by MIR spectra depends on the coupling of the dust grains to the gas in the inner disc, as the dust grains that deliver the volatiles provide additional continuum opacity that can hide the vapour; and this is particularly true for water. As such, there is no one-to-one relationship between column density and vapour mass. Thus, unless the amount of dust and the extinction it produces can be directly quantified, absolute column densities should be interpreted with caution, as elevated values could reflect delivered material or a loss of dust (Fig. 7). Tentatively, the relatively low column densities observed with MIRI-MRS could reflect a high retention of dust in the inner disc (Scenario 1), but they could also indicate that the observed spectra are dominated by a thinner layer of hotter material at the disc surface.

3. On the other hand, a ratio of column densities – and in particular the CO₂ /H₂O ratio $N_{{\text{CO}}_2}/N_{{\text{H}}_2\text{O}}$ – more closely reflects the underlying chemistry of the disc and can thus be used to trace the delivery of molecules to the inner disc. We predict this ratio to be enhanced on Myr timescales as a signature of this sequential drift and advection, in agreement with some observed CO₂-dominated spectra (Fig. 8, grey lines), but sensitivity to weak spectral features that probe deep into the disc is likely needed to be able to see it.

4. Dust traps can significantly alter the molecular abundances of the inner disc by interrupting the radial drift of dust and the supply of molecules. As a result, the peak enhancements are reduced and also occur somewhat earlier. As the H₂O vapour is accreted on a shorter timescale than the CO₂, it is more quickly affected by the lack of resupply, which increases the CO₂/H₂O ratio (Fig. 8, blue lines).

5. Assuming they exist from very early in the disc’s evolution, closer-in traps block more dust, creating an initially positive dependence of molecular abundances on the gap location. This switches to a negative dependence, driven by the leakiness of the traps after several viscous timescales (at the snow line). The lag of the CO₂ evolution behind that of H₂O means that their ratio also has a positive dependence on gap location at t ≲ 2 Myr; however, at later times, the two leak from the trap at the same rate (Fig. 9, bottom panels), thus erasing any dependence.

6. Adding a delay to the formation of a trap most strongly affects close-in traps, where the growth and drift timescales are short. Formation on timescales of up to 1 Myr still prevents the long-term resupply of H₂O ice on Myr timescales, while letting more CO₂ ice through initially, thus raising the CO₂/H₂O ratio. This reduces the strength of the dependence of this ratio on the trap location (Fig. 10).

Therefore, we suggest that if carefully measured, $N_{{\text{CO}}_2}/N_{{\text{H}}_2\text{O}}$ retrieved from MIR spectra of protoplanetary discs can be used as a tracer of radial drift. Moreover, the presence of a positive relationship across a population of discs between $N_{{\text{CO}}_2}/N_{{\text{H}}_2\text{O}}$ and the location (or at least presence) of a gap could be considered as evidence of the influence of dust traps on this drift, with its strength informing us about the timescale on which the traps form. We suggest that, based on current observations of or constraints on their substructures, the CO₂-dominated discs GW Lup and CX Tau may qualitatively agree with the expected trend.

Data availability

The version of the DiscEvolution code used for our models can be found at https://github.com/AndrewSellek/DiscEvolution/releases/tag/v2.0.0.

Appendix C can be found on zenodo at https://doi.org/10.5281/zenodo.14269071.

Acknowledgements

We thank the referee for their useful input on the comparisons between our models and different spectral features. The authors thank members of the MINDS team – including T. Henning, I. Kamp, D. Gasman, S. Grant, N. Kurtovic, and M. Temmink – for stimulating discussions and also thank M. Temmink for advanced information about slab model fits for DR Tau. We are also grateful to R. Booth for developing the DiscEvolution code and for discussions about implementing gaps therein. A.D.S., M.V., and E.v.D. acknowledge support from the ERC grant 101019751 MOLDISK. E.v.D. also acknowledges support from the Danish National Research Foundation through the Center of Excellence “InterCat” (DNRF150) and grant TOP-1 614.001.751 from the Dutch Research Council (NWO).

Appendix A Transport vs chemical processing

Here we set out to validate the assumption that for the main species considered in this work we can ignore chemical reactions. Figure A.1 shows the rate on which the vapour mass M_vap,X is changing, defined as $ℛ=\frac{\text{d}}{ \text{d}t}\ln \left(M_{\text{vap },\text{X}}\right)$, for both H₂O and CO₂ in a smooth disc for each scenario. Scenarios 1 and 3 present very similar evolution in this space as expected from the fact that their differences should mainly affect the retention of dust inside the H₂O snow line, and therefore only quantities relating to the amount of observable molecules.

The initial evolution is characterised by the growth and drift of dust and delivery of volatiles. As drift timescales for fragmentation-limited dust scale with radius t_drift ∝ R^1/2, and growth timescales scale roughly as $t_{\text{grow }}~\frac{10}{ϵ\text{Ω}}\propto R^{3/2}$ (assuming roughly 10 e-foldings are needed to grow the dust from its initial size to when it starts to decouple), for the bulk of disc material at large radii, its ability to supply of volatiles is limited mainly by the growth timescale of the dust. The mass of dust that has grown to pebble size at time t is then that contained within the pebble production front R_grow(t) ∝ t^2/3 and M(< R) ~ R/R_C for R < R_C; the rate of change of the vapour mass is thus $~\frac{2}{3t}$. Rates of 10⁻⁵ − 10⁻⁴ yr⁻¹ are seen on timescales of 10⁴ − 10⁵ yr, consistent with this estimate.

Once the first generation of pebbles has passed through the disc, the pebble flux drops and the disc switches to draining the molecules through viscous accretion. The rate, here is simply the inverse of the timescale given by Eq. 19. As a consequence, H₂O is removed at a rate of ~ 10⁻⁵ yr⁻¹, and CO₂ at a rate of ~ 2 × 10⁻⁶ yr⁻¹. These are in good agreement with the most negative rates seen in Fig. A.1. Scenario 3 reaches rates approximately 10 times smaller, consistent with its 10 times lower value of α.

Included for comparison is the photodissociation rate for CO₂ due to the UV field produced by cosmic rays, which is typically around 10⁻¹⁴ s⁻¹ (Bosman et al. 2018). This is typically faster than the rate at which cosmic-ray-ionised He can destroy molecules. The delivery rates in the first ~ 10⁵ yr are comfortably above these expected ≲ 3 × 10⁻⁷ yr⁻¹ (photo-)chemical destruction rates of molecules and therefore we are justified to ignore (photo-)chemical reactions at the midplane during this period. Similarly, the removal of molecules from the inner disc also happens much more quickly by gas accretion than by chemical destruction, at least for the first ~ 5 Myr. The (photo- )chemical reactions can therefore only become relevant for H₂O and CO₂ on the longest timescales of disc evolution, when the net timescale resulting from the balance of removal by accretion and delivery by drift tends towards smaller values because the molecular abundances are low enough to be efficiently replaced by ices sublimating from the small grains entrained with the gas. However, Scenario 3 is likely somewhat more susceptible to the removal of CO₂ later in evolution happening via chemical reactions rather than accretion due to its low viscosity.

Overall, this suggests that indeed for most of a disc’s evolution (except at low viscosities), the transport processes are faster than typical chemical processing timescales for the main molecules of interest in this work - at least in the vertically integrated sense that a 1D model such as ours explores - and thus we are justified in neglecting (photo-)chemical reactions. This is consistent with the results of Booth & Ilee (2019), who saw that transport of molecules to the inner disc wiped out the signatures of chemical transformation for H₂O and CO₂.

We note that the impact of faster chemical processing at the disc surface (where molecules are exposed to higher UV fluxes) will depend on how the resulting reaction rates compare to the vertical mixing timescale. Assuming a turbulent model for the vertical diffusion with α_z = 10⁻³, the inverse timescale will be ~ 0.01 yr⁻¹ at the ~ 0.5 au effective emitting radius of the emission. We refer to Semenov & Wiebe (2011) for a more detailed discussion of the comparison between vertical mixing timescales and the rates of different types of reaction.

Appendix B Trends with gap location for scenarios 2 and 3

In Scenario 3 (Fig. B.4), the trends in the vapour masses strongly resemble that of Scenario 1 (Fig. 9), except for the fact the negative trend gets even steeper. In Scenario 2 (Fig. B.3) however, the evolution of the disc gas is much slower and hence the trends do not undergo a reversal during the 10 Myr modelled lifetime. This slower evolution also means that the trends for H₂O are much closer to those predicted by Eq. 23.

We also note that Scenario 3 shows a steeper negative relationship than the other two. This is because in this low- turbulence scenario, fragmentation due to differential drift limits the grain sizes (apart from at the very centre of the trap). In this case, F_leak is modified to $$F_{\text{leak }}\propto \frac{R_{\text{gap }}}{w_{\text{gap }}^2}\exp \left(-\frac{R_{\text{gap }}}{R_{\text{C}}}\right).$$(B.1)

Since 8 predicts that $\frac{\alpha }{St}\propto c_{\text{S}}^2/v_{\text{K}}\propto R_{\text{gap }}^0$, then w_dust ∝ w_gas ∝ R^5/4, resulting in a steeper F_leak ∝ R^−3/2.

Moreover in this case, as one barrier to dust growth is removed, the dust can grow somewhat larger and become more concentrated, making the diffusive flux even stronger and the traps generally leakier. This is somewhat counter-intuitive but results from the assumption of a fixed f_small in the two-population model; a more realistic model of coagulation and fragmentation would likely suggest a smaller f_small and a net reduction in the leaking flux. Therefore, while the predicted steeper trend is probably correct, the relative values with respect to Scenario 1 may not be. Indeed, an alternate fate for the dust trapped outside the gaps that—which would reduce the flux of dust leaking into the inner disc— is instead to become converted to planetesimals due to the large concentration of solids. Zagaria et al. (2023) proposed that the outer dust ring of HD HD 163296 may be undergoing planetesimal formation by the Streaming Instability (SI) and planetesimal formation has also been proposed to regulate the optical depth of dust rings, explaining the typical values observed at (sub)millimetre wavelengths (Stammler et al. 2019). However other works have argued that pressure bumps may not be suitable sites for SI to operate to turn millimetre-sized grains into planetesimals (Carrera et al. 2021; Carrera & Simon 2022). Kalyaan et al. (2023) included planetesimal formation in their models, generally finding it only to be triggered within the traps unless the fragmentation velocities were very high, in which case some planetesimal formation also occurs at the snow line. As they focused on the earlier stages where the H₂O was enhanced in the inner disc, Kalyaan et al. (2023) argued that these solids were trapped anyway, and hence on these timescales there were negligible differences between models with and without planetesimals. Regardless, so long as the trap is outside both the CO₂ and the H₂O snow line, both ices should be equally incorporated into planetesimals and thus while the late time abundances may fall, this should not affect the CO₂/H₂ O ratio.

Unlike Scenario 1, the lack of a traffic jam inside the H₂O snow line in both Scenarios 2 & 3 gives us more variation in H₂O column densities. On timescales of ≳ 2 Myr, we do see clear positive dependences of the column densities on R_gap in Scenario 2; however, as previously discussed, the column densities of both CO₂ and H₂O are very high compared to observed values. A further consequence of such a scenario would be that any of the lines commonly used to measure H₂O emission (Gasman et al. 2023, 2025) would likely saturate and provide no sensitivity to the underlying abundances.

Fig. A.1 Evolution of the rate at which the vapour masses of H2O (left) and CO2 (right) are changing in a smooth disc model in Scenarios 1 (purple), 2 (gold), and 3 (red) as a result of delivery by drifting dust and viscous advection onto the star. A typical rate for cosmic ray (CR) photodissociation of CO2, 10−14 s−1, is shown as the dashed lines as a reference for the rates at which chemical reactions occur.

In Scenario 3, we actually see both the positive and negative dependences reflected in the column densities (at early and late times respectively) for both CO₂ and H₂O. The ability to see the late-time trends suggests that the dust that leaks through the trap is dynamically decoupling from the gas, because the leakier traps promote a higher inner disc dust abundance and therefore more frequent collisions leading to grain regrowth in the inner disc. In the radial-drift–limited regime, the dust surface density becomes dependent on $F_{\text{leak }}^{1/2}$ (c.f. Birnstiel et al. 2012, Eq. 39), whereas the vapour surface density - or the dust surface density when the dust remains coupled - scales as F_leak (c.f. Birnstiel et al. 2012, Eq. 41). Therefore, in Scenario 3, applying the approximation of Eq. 22, the visible column density scales as $F_{\text{leak }}^{1/2}$ thus introducing a negative dependence on R_gap. For the closest-in gaps, the dust growth is so effective that it reaches the fragmentation barrier, whence it again depends linearly on F_leak and the contributions of the vapour and dust to the visible column density again cancel. Therefore, a negative relationship between the column density ratio and the location of a dust trap observed across a population of disc may indicate that dust is leaking through their trap and efficiently regrowing to drifting sizes.

Fig. B.1 Same as Fig. 8 but for Scenario 2, which has no traffic jam inside the water snow line and a slow gas accretion.

Fig. B.2 Same as Fig. 8 but for Scenario 3, which has no traffic jam inside the water snow line and a faster gas accretion.

Fig. B.3 Same as Fig. 9 but for Scenario 2, which has no traffic jam inside the water snow line and a slower gas accretion.

Fig. B.4 Same as Fig. 9 but for Scenario 3, which has no traffic jam inside the water snow line and a fast gas accretion.

References

All Tables

All Figures

Fig. 1 Cartoon depicting the three scenarios explored in this work. The thickness of the yellow arrows connecting the disc and star represents the magnitude of α controlling the gas evolution and accretion. The turbulence strength is indicated by the thickness and size of the circular arrow in the bottom right. The speckled region represents the dust with the brown colour indicating dry dust and the light blue colour indicating icy dust and the finer-grained pattern indicating where grains are more fragile. The vertical size of this region is inflated for illustrative purposes in locations where is a traffic jam and grains pile up.

In the text

Fig. 2 Example surface density profiles of gas (solid, initial value only) and (small+large) dust (dashed, at different times) for each of the three scenarios, illustrating the changing dust surface density distribution due to radial drift. The dip in the gas surface density signifies the gap, in this case at Rgap = 30 au and after 0.1 Myr (red), dust has begun to accumulate in the trap at its outer edge. The traffic jam that occurs in Scenario 1 can be seen in the steep rise in the dust surface density after 0.1 Myr inside the water snow line at ≈0.4 au. The other small steps correspond to changes in the solid fraction at snow lines and the smaller bumps are the accumulation of ices at snow lines due to the cold finger effect.

In the text

Fig. 3 Evolution of the pebble component of the disc. Right: total dust mass (including that which has crossed the inner disc boundary and accreted onto the star) inside the initial locations of the H2O (solid) and CO2 (dashed) snow lines (right). Left: Pebble flux across the H2O snow line. Included for reference are the values from Table 6 of Romero-Mirza et al. (2024a) [RM24], which are estimated from the total number of H2O molecules <400 K, split into compact discs (small, filled, markers) and extended discs (larger, ringed, markers).

In the text

	Fig. 4 Radial distribution of the H2O (blue), CO2 (bright green), and CH4 (olive green) after 1 kyr of evolution (left) and 0.1 Myr of evolution in a smooth disc in Scenario 1. Solid lines indicate molecules in the solid (ice) phase, and dashed lines molecules in the gas (vapour) phase.
In the text

	Fig. 5 Evolution of the vapour mass of H2O (blue), CO2 (bright green), and CH4 (olive green) normalised to their initial values for a smooth disc in Scenario 1. The vertical dotted lines mark the times that each molecule reaches its peak vapour mass and the rightward arrows indicate the timescale for the peak to then decrease (Eq. (19)).
In the text

Fig. 6 Example of the construction of a synthetic spectrum and the retrieval of parameters using 0D LTE slab model fits. The top panels show the temperature profile (left) and the column density profiles after 104 yr (right) for small dust (black) and total (solid) and visible (dashed) H2O (blue) and CO2 (green). The light, vertical, dotted lines indicate the snow lines. The middle panel shows the synthetic spectra of each molecule (assuming a distance of 140 pc) while the bottom panels show the goodness of fit for a grid of slab models with varying temperature and column density. The white contours show the radius (in AU) of a circle whose area maximises the goodness of fit for each combination of temperature and column density. The overall best fit parameters are shown on the top panels as the crosses, with the vertical error bars representing the spacing of the grid of slab models.

In the text

	Fig. 7 Relationships between the vapour mass and the column density measured using slab fits for H2O (left) and CO2 (right) for a smooth disc model in Scenarios 1 (purple), 2 (gold), and 3 (red). The curves shown are traversed anticlockwise.
In the text

Fig. 8 Evolution of the vapour mass (left) and column density (right) for CO2 (top row), H2O (middle row), and their ratio (bottom row) for discs with traffic jams inside the water snow line (Scenario 1). The grey line represents a smooth disc with no traps, while the blues lines, in order of increasingly dark shade, represent traps with increasingly close in traps. Overplotted at their approximate age are the results of slab model fits for a few discs from the literature: CX Tau (Vlasblom et al. 2025), GW Lup (Grant et al. 2023), SY Cha (Schwarz et al. 2024), Sz 98 (Gasman et al. 2023), DR Tau (Temmink et al. 2024b,a), AS 209 (Romero-Mirza et al. 2024b) and DoAr 33 (Colmenares et al. 2024), ordered by increasing stellar mass (more yellow colours). The open squares indicate an alternative estimate obtained from the 13CO2 isotopologue. The grey bar in the top left shows the ± log10(2.5) uncertainty depending on the basis used for the initial ice abundances (Sect. 2.3.2).

In the text

Fig. 9 Trends of CO2 (top) and H2O (middle) vapour masses (left) and column densities (right) and their ratio (bottom) with Rgap at different stages of evolution for discs with traffic jams inside the water snow line (Scenario 1). The initial condition is in grey, while increasingly dark shades of green correspond to increasingly late times. The model with Rgap = 1 au is plotted separately at the dots as its gap is located very close to the H2O snow line and thus inside the CO2 snow line. A smooth disc model is included for comparison as the triangles. Equations (23) (dashed-dotted) and (11) (dotted) are plotted as guides to the eye the trends for approximate trends at early and late times.

In the text

	Fig. 10 Dependence of the CO2/H2O column density ratio on gap location at different snapshots in the disc evolution. The panels show models in which the traps are fully developed at successively later times.
In the text

Fig. A.1 Evolution of the rate at which the vapour masses of H2O (left) and CO2 (right) are changing in a smooth disc model in Scenarios 1 (purple), 2 (gold), and 3 (red) as a result of delivery by drifting dust and viscous advection onto the star. A typical rate for cosmic ray (CR) photodissociation of CO2, 10−14 s−1, is shown as the dashed lines as a reference for the rates at which chemical reactions occur.

In the text

	Fig. B.1 Same as Fig. 8 but for Scenario 2, which has no traffic jam inside the water snow line and a slow gas accretion.
In the text

	Fig. B.2 Same as Fig. 8 but for Scenario 3, which has no traffic jam inside the water snow line and a faster gas accretion.
In the text

	Fig. B.3 Same as Fig. 9 but for Scenario 2, which has no traffic jam inside the water snow line and a slower gas accretion.
In the text

	Fig. B.4 Same as Fig. 9 but for Scenario 3, which has no traffic jam inside the water snow line and a fast gas accretion.
In the text