© The Authors 2025

1 Introduction

The inner regions (<10 au) of planet-forming discs are active sites of planet formation (Morbidelli et al. 2012; Dawson & Johnson 2018). The chemical composition of these forming planets is set by the elemental and molecular composition of the nascent disc. One of the key ingredients for habitable worlds is H₂O and, therefore, studying the available H₂O vapour reservoirs in discs is of great importance.

Based on H₂O vapour observations with the Spitzer Space Telescope and spatially resolved continuum images with the Atacama Large Millimeter/submillimeter Array, Banzatti et al. (2020) proposed scenarios for the expected H₂O reservoirs given the sizes of the dust discs. Discs that are observed to be compact in their millimetre emission (R_dust < 60 au) are thought to have efficient radial drift, leading to enriched H₂O reservoirs (N ∼ 10¹⁸–10²¹ cm⁻²) in the innermost regions (<few au) of these discs. The larger, structured discs (60< R_dust < 300 au) are expected to not show these enriched H₂O abundances (≤ 10^17.5 cm⁻²), since sub-structures (i.e. gaps and rings) could trap the icy dust grains in the outer regions, halting them from reaching the inner disc and from crossing the H₂O snowline. Furthermore, for discs with large cavities, the situation is expected to be even more extreme, with the inner discs expected to be depleted of H₂O.

Recent modelling studies have explored the notion of both unimpeded enhancement due to radial drift and the opposite case of sub-structures halting the inward drift of pebbles. In particular, Kalyaan et al. (2023) investigated the influence of gaps on the gaseous enhancement in the inner disc. They noted that the enhancement, at least for the H₂O vapour reservoir, is of limited duration (up to a few million years) before the material gets accreted onto the host star. It is only when the gaps do not block the dust entirely, but some of the grains are still able to pass through (Pinilla et al. 2024; Mah et al. 2024) that the lifetime of this enhanced reservoir be prolonged. However, the inward drift of the dust particles not only enhances the gaseous reservoir, but it also increases the opacity of the dust itself. Recent models by Sellek et al. (2025) suggest that the increased opacity of the dust elevates the τ=1 layer of the continuum, which could hide a larger gas reservoir that lies deeper in the disc. The observable column densities may thus not reflect the enhanced reservoirs. Further works by Houge et al. (2025) support this idea, while also investigating the competition between photodis-sociation and vertical mixing on the available reservoir. Finally, Kalyaan et al. (2023), Sellek et al. (2025), and Easterwood et al. (2024) noted the importance of the gap locations and the time at which they form for the drifting pebbles. Gaps located at smaller radial distances, and which formed early on, are more effective at limiting the gaseous enhancement, as more dust grains are blocked outside the snowline. However, they may leak more for a given gap depth or grain size. For gaps located at larger radial distances, the enhancement is stronger, as a smaller icy dust reservoir is blocked. Furthermore, the way the gap opens may also be of importance. Lienert et al. (2024) have shown that gaps opened due to internal photoevaporation are able to block both the pebble and gas flow, strongly influencing the inner disc composition. In contrast, the models from Greenwood et al. (2019) have shown that in the absence of gaps, the dust opacity decreases as the disc evolves with time and this leads to an increase in the H₂O flux. Current observations of mid-infrared spectra of large and structured discs with the James Webb Space Telescope (JWST), including those with large cavities, already show that inner regions of such discs are not completely devoid of H₂O (Perotti et al. 2023; Schwarz et al. 2024). Furthermore, in some cases, they may have strong emission from the cold reservoir (<400 K; Gasman et al. 2023, 2025).

The James Webb Space Telescope (Rigby et al. 2023) provides the best opportunity to fully explore the available H₂O reservoirs in the inner regions of planet-forming discs, using the increased sensitivity and resolution of the Medium Resolution Spectrometer (MRS; Wells et al. 2015; Argyriou et al. 2023) of the Mid-InfraRed Instrument (MIRI; Wright et al. 2015; Rieke et al. 2015; Wright et al. 2023) with respect to the Spitzer Space Telescope. Since its launch, multiple works have studied the available reservoirs, including the analysis of H₂O emission across the entire JWST-MIRI wavelength range (Gasman et al. 2023) using 0D local thermal equilibrium (LTE) slab models. This analysis directly proved the existence of an expected radial temperature gradient, where the longer wavelengths probe larger radii and, thus, colder temperatures (Blevins et al. 2016; Banzatti et al. 2017, 2023b, 2025). Banzatti et al. (2023a) identified the emergence of a cold H₂O reservoir (T < 240 K), which is expected to be the effect of radial dust drift. Subsequently, multiple slab models were fitted to the pure rotational H₂O spectra (> 10 µm) identifying the available reservoirs (Pontoppidan et al. 2024; Temmink et al. 2024a). A third cold component, in addition to hot (> 800 K) and warm (400 < T < 800 K components, was found to be necessary to describe the same cold H₂O reservoir of T < 240 K. The three-component analysis of Temmink et al. (2024a) also provided another confirmation of the radial temperature gradient, where the multiple-component analysis of DR Tau could be approximated with a temperature profile of T(R) ∼ 500 (R/1au)^−0.5 K. More recently, efforts have been made to increase the complexity by moving away from using multiple components and to describe the temperature and column density profiles by parametric functions (Kaeufer et al. 2024; Romero-Mirza et al. 2024). These profiles, including simple and exponentially tapered power laws, show that the pure rotational H₂O transitions can be explained by such parametric models; however, this has only been tested for a limited sample of discs with different disc characteristics, such as their radial sizes and structures. We also note the effort by Woitke et al. (2024), who used a full 2D thermochemical code to model, for the first time, the observed molecular emission in the outbursting source of EX Lup.

In this work, we aim to analyse the pure rotational H₂O emission in 8 millimetre-compact discs (dust radii of R_dust,95% ∼ 25–60 au; Long et al. 2019; Facchini et al. 2019), using both the parametric and multiple component techniques highlighted above. Our criterion for a dust disc to be compact follows the classification of Banzatti et al. (2020), namely, assuming compact discs have dust radii of ≲60 au. All sources have been observed with JWST/MIRI-MRS as part of the JWST Guaranteed Time Observations Program MIRI mid-INfrared Disk Survey (MINDS, PID:1282, PI: T. Henning; Kamp et al. 2023; Henning et al. 2024). Throughout our analysis, we aim to investigate the strengths of the different H₂O reservoirs, indicated throughout the paper as hot (>800 K), intermediate (400–800 K), and cold (<400 K, see also Banzatti et al. 2023a). We use this to explore the H₂O reservoirs and column density profiles. In addition, we also infer whether these compact discs show an enhancement in the cold H₂O reservoir, as suggested by Banzatti et al. (2020), or whether the situation is more complex than this scenario.

The paper is structured as follows. Section 2 describes the sample and the observations, while Section 3 contains a description of the used methodology. The results are represented in Section 4 and interpreted in Section 5. Aside from the interpretation, Section 5 also discusses the role of sub-structures in setting the H₂O reservoirs and offers a comparison with existing models. Finally, Section 6 contains our conclusions and a short summary.

2 Sample and observations

2.1 Sample

Our sample of compact planet-forming discs consists of eight sources: BP Tau, CX Tau, CY Tau, DN Tau, DR Tau, FT Tau, RNO 90, and XX Cha, for which stellar properties are summarised in Table 1. Dedicated papers, exploring the observed molecular emission, exist for two of the sources: CX Tau (Vlasblom et al. 2025) and DR Tau (Temmink et al. 2024b,a).

Although the discs in our sample are known to be rather compact in the millimetre continuum emission (R_dust ≲ 60 au), there is limited information available about their gas disc sizes (R_gas). Trapman et al. (2019) proposed that a gas-to-dust size ratio of R_gas/R_dust > 4 is a clear sign of radial drift. For our sources, only one disc has literature values for both R_dust and R_gas: CX Tau (Facchini et al. 2019). These values suggest a ratio of R_gas/R_dust > 5, hinting at strong radial drift being the potential cause of this disc’s compactness. For the other discs, the gas radius has not been determined and, in some cases, high-resolution observations of the continuum emission do not (yet) exist. Consequently, the gas-to-dust size ratio is still to be determined in a consistent manner and it is unclear whether the compactness of these discs can be attributed to an evolution dominated by radial drift or simply to their initial conditions. For these reasons, we have chosen our sample to follow the criteria of Banzatti et al. (2020) and we considered discs with millimetre continuum emission radii of R_dust ≲ 60 au as compact.

As can be seen from Table 1, our sample is rather homogeneous in terms of masses, luminosities, and inclinations. There are a few outliers: for instance, DR Tau and RNO 90 have higher stellar masses compared to the others, while DR Tau also has the lowest inclination (i ∼ 5.4^°, Long et al. 2019). CX Tau, on the other hand, has the highest inclination (∼ 55.1^°). Furthermore, as six out of eight discs are located in the Taurus star-forming region, we can expect them to be of similar ages (∼ 1–3 Myr, Krolikowski et al. 2021; Luhman 2023). In contrast, RNO 90 is located in Ophiuchus and may, therefore, be on the younger side. These small differences suggest that our sample is overall rather homogeneous and differences must therefore be due to either their initial conditions or differences in their (potentially drift-dominated) evolution. A larger sample will be analysed in Temmink et al. (in prep.), which may provide more insights into the differences related to the stellar properties.

2.2 Observations and data reduction

The JWST-MIRI/MRS details (date of observation and integration time) are also included in Table 1. All MRS spectra have been taken in the FASTR1 readout mode with a four-point dither pattern using all three grating settings (A, B, and C). All data have been reduced using a standard pipeline reduction (version 1.16.1; Bushouse et al. 2024) and pmap 1315. The spectra have been extracted through aperture photometry, where the aperture has a size of 2× the full width at half maximum (FWHM). Additionally, we have corrected for residual fringes using the implementation of the default pipeline, as most of the observations (except XX Cha) were taken without target acquisition. The resulting spectra were continuum subtracted using an updated version of the method by Temmink et al. (2024b), which makes use of the ‘iterative reweighted spline quantile regression’ method included in the PYBASELINES package (Erb 2022). Before estimating the continuum, downward spikes were masked using the same method as in Temmink et al. (2024b), but now they were masked per MIRI sub-band as opposed to over the whole spectrum. We used a quantile regression value of 0.1 and placed the knots of the cubic splines every 75 data points, except for Channel 4 and the silicate feature (∼ 8.25–11.25 µm), where the knots were placed every 25 points to better estimate the varying continuum. For CX Tau and DN Tau, a knot spacing of 25 data points was used for the entire spectrum, due to their strongly varying continuum. Given the large uncertainty of the observations at the location of the silicate feature, the quantile regression value was set to 0.5 for all discs, which places the baseline through the median of the observations. Additionally, we changed the continuum estimate for the Q-branches of CO₂, HCN, and C₂H₂. The typical fit could overestimate the continuum of these Q-branches, effectively taking away molecular flux in the subtraction. To avoid this over-subtraction, we used a cubic interpolation of the estimated baseline just before and after the Q-branch to better fit the continuum level underneath the Q-branch. In particular, we masked the 13.45–14.20 µm wavelength region that captures the HCN and C₂H₂ Q-branches, and the 14.88–15.01 µm region for that of CO₂. The reduced spectra and continuum estimates are displayed in Figure 1.

Figure 1 shows that the shape of the spectra overall looks similar and in some cases nearly indistinguishable. All sources have a silicate feature at ∼ 10 µm and contain emission features from a variety of molecular species. The only outlier is CY Tau, which has stronger emission at the shortest wavelengths. This is likely due to an inner (puffy) rim self-shadowing the outer disc and, subsequently, lowering the mid- and far-infrared flux (see, for example, Dullemond & Dominik 2004; Woitke et al. 2019). We leave the analysis of the dust continua to a future work.

3 Methodology

To analyse the pure rotational H₂O emission in the JWST-MIRI/MRS spectra of our sample of compact discs, we used slab models under the assumption of local thermal equilibrium (see Grant et al. 2023; Tabone et al. 2023 for more details on the slab model generation). Following Banzatti et al. (2025), we included mutual line shielding in the H₂O slab models to account for the shielding of ortho- and para-line pairs. Additionally, we took the mean spectral resolution, calculated using the results of Pontoppidan et al. (2024) and the updates by Banzatti et al. (2025), of each MIRI sub-band to properly sample the wave-length grid of the slabs over the entire spectrum. We used two different line widths for our models: a constant value of ΔV=4.71 km s⁻¹, which is the line width of H₂ at a temperature of 700 K (Salyk et al. 2011), as well as a similar approach to Romero-Mirza et al. (2024), where the line width is taken as the sum in quadrature (hereafter, simply called the quadrature line width) of the turbulent line width (ΔV_turb fixed to 1.0 km s⁻¹) and the thermal line width ($\text{Δ}{V}_{\text{therm }}=\sqrt{2kT/m}$, with T the temperature of the slab model and m the mass of an H₂O molecule). The two line widths were used to make a comparison between the practices applied in the recent literature. While earlier works (e.g. Tabone et al. 2023; Grant et al. 2023; Gasman et al. 2023) have used a constant value of 4.71 km s⁻¹, Romero-Mirza et al. (2024) recently used the sum in quadrature. To investigate which of these approaches offers a better fit to our spectra, we tested both options in our fits and make a comparison (see Section 5.4). The Python package SPECTRES (Carnall 2017) was subsequently used to resample the slab models onto the wavelength grid of MIRI. However, before the resampling, we applied a wavelength shift to the slab models (see also Pontoppidan et al. 2024), based on the heliocentric velocity of the sources (see Table 1). Figure A.1, using the spectrum of DR Tau as an example, shows that these wavelength shifts, although subtle, are required to properly recreate the observed line profiles and will subsequently improve the fit results.

We fit our spectra using two different approaches: similarly to Temmink et al. (2024a), we used a multi-component analysis where three or two distinct components (with decreasing temperatures and increasing emitting areas) were fitted to the entire rotational spectrum. Furthermore, we use and extended the methods implemented by Romero-Mirza et al. (2024), who used (exponentially tapered) power law profiles for the temperature and column density. We denoted these models as multi-component and parametric models; in the latter case, we explored more options than (exponentially tapered) power laws. By using both these methods, we have been able explore the strengths of the different H₂O reservoirs and investigate the robustness of what kind of profile for the column density is preferred in the inner region. By using both techniques, we can also gain more insights into the changing excitation conditions with radius in the inner regions of these discs.

Fig. 1 Normalised-to-peak-flux spectra of our sample of millimetre-compact discs. The red line indicates the estimated continuum. The indicated types and their meaning are discussed in Section 5.1. In blue, we highlight the 23.72–24.03 µm wavelength region, where two transitions are located that are most important for obtaining information about the cold H2O reservoir.

Fig. 2 Schematic highlighting the different parametric models used for the column density profiles. Both the horizontal (radius) and vertical (column density) axes are taken in log10-space.

3.1 Rotational H₂O spectrum: Parametric analysis

In an attempt to parametrically describe the observed H₂O emission, we followed the approach from Romero-Mirza et al. (2024), whereby 50 slab models are used to sample profiles in both temperature and column density as a function of the emitting radius. As opposed to Romero-Mirza et al. (2024), we chose to first determine the radius as a function of temperature, $R(T)=R_{\text{in }}{\left(\frac{T}{1500\text{K}}\right)}^{-1/q}.$(1)

This relation ensures that the inner radius (R_in (in au), as listed in Table 1) of our slab models is located at the estimated dust sublimation radius at 1500 K. Using this equation, we sampled the annular emitting regions of our 50 slab models for temperatures between 1500 K (the dust sublimation radius; Barvainis 1987) and 150 K (the condensation temperature of H₂O; Collings et al. 2004) in log₁₀-space. The dust sublimation radius was calculated following the approach of Dullemond et al. (2001), $R_{\text{in}}\sim 0.07\sqrt{\left(L_{\ast }/L_{\odot }\right)}$. Both the calculated inner radii and adopted stellar luminosities can be found in Table 1. By simply rewriting Equation (1), we obtain a relation for the temperature as a function of the radius, $T(R)=1500\hspace{0.17em}\text{K}{\left(\frac{R}{R_{\text{in}}}\right)}^{-q}.$(2)

For the column density, we used the same profiles as Romero-Mirza et al. (2024): a power law (profile I, Equation (3)) and an exponentially tapered power law (profile II, Equation (4)). In addition, we also tried different profiles. Figure 2 provides a schematic overview of the different profiles and their potential behaviours. In particular, we assumed a constant number of molecules (𝒩_mol) with radius, whereby each annulus contains the same number of molecules, which automatically resembles a power law with a power of –2 (the emitting area, A, is proportional to the square of the radius). We then allowed this number of molecules to experience a jump; namely, the number of molecules would get enhanced (or depleted) by a factor, F_scale, at temperatures below a given jump temperature, T_jump (see Equation (5)). We tried profiles with a fixed jump temperature, at 400 K (profile III) and we kept T_jump as a free parameter (profile IV). The value of 400 K approximately resembles the boundary temperature between the sublimated (colder temperatures) and gas-phase-formed (hotter temperatures) H₂O vapour reservoirs (see also Romero-Mirza et al. 2024). The jump temperature can move to higher temperatures, indicating a potential depletion of the hot H₂O reservoir, or to colder temperatures, representing a potential enhancement of the cold H₂O reservoirs. A parabola provides a smoother description of this behaviour: the jump moving to higher temperatures is captured by an upward parabola (the parabola has a maximum), while a downward parabola (the parabola has a minimum) captures the behaviour of the jump moving to lower temperatures (see also Figure 2). Therefore, we fitted a parabola in log space (profile V, Equation (6)). The equations describing the column density profiles of each of these cases are: $N_{\text{I}}(R)=N_0{\left(\frac{R}{1\hspace{0.17em}\text{au}}\right)}^{-p},$(3) $N_{\text{II}}(R)=N_0{\left(\frac{R}{1\hspace{0.17em}\text{au}}\right)}^{-p}\exp \left[-{\left(\frac{R}{R_{\text{c}}}\right)}^{\varphi }\right],$(4) $N_{\text{III,IV }}(R)\propto \{\begin{array}{cc}\frac{{𝒩}_{\text{mol }}}{A}& \text{if} T>T_{\text{jump }},\\ \frac{{𝒩}_{\text{mol }}\times F_{\text{scale }}}{A}& \text{if} T\le T_{\text{jump }},\end{array}$(5) $N_{\text{V}}(R)=N_0{\left(\frac{R}{1\hspace{0.17em}\text{au}}\right)}^{\alpha +\beta \ln (R/1\text{au})}.$(6)

The emitting area of each slab model is accounted for in the exact same manner as reported in Romero-Mirza et al. (2024). The area is multiplied by the cosine of the disc’s inclination (i, reported in Table 1), $A=\pi \left(R_2^2-R_1^2\right)\times \cos \hspace{0.17em}(i)$, where R₁ and R₂ are, respectively, the inner and outer radii of each annulus (see also Equation (1)).

We fit the profiles, sampled by the 50 slab models, using a Markov chain Monte Carlo (MCMC) approach, implemented through the EMCEE python-package (Foreman-Mackey et al. 2013). From the posterior distribution of the MCMC, we take the median values as the best-fit values for the fitted parameter. Additionally, the 18^th and 84^th percentiles are taken as the lower and upper 1σ uncertainties. We use a reduced-χ² as the likelihood function of the MCMC, in order to compare the full model (the sum of the 50 slab models) with the observations. The likelihood has, therefore, takes the following form: ${\chi }_{\text{red}}^2=\frac{1}{N_{\text{data}}}\underset{i}{{{\displaystyle \sum }}^{\text{}}}\frac{w_i{\left(F_{\text{obs},i}-F_{\text{model},i}\right)}^2}{{\sigma }_i^2},$(7)

where F_obs is the observed spectrum and N_data is the number of data points within the used fit regions. Then, σ denotes the noise of the observations, determined in a similar manner as in Temmink et al. (2024b), where the James Webb Exposure Time Calculator¹ has been used to obtain the continuum signal-to-noise (S/N) at the midpoint of each sub-band given the observational setup. The S/N values are given in Table B.1. As the noise is significantly higher in Channel 4 and many important H₂O lines can be found at these wavelengths, we slightly offset the large noise by introducing weights in the χ²-formula, w_i. All fit regions beyond 20 µm have been given a higher weight of 5, while select regions with important transitions probing the cold reservoir have been given a weight of 10 or 15 (for the strongest transitions). While the values are arbitrary, these weights ensure that we are able to properly fit the full rotational spectrum, without favouring the hot reservoir, as the transitions at smaller wavelengths have lower flux uncertainties. Fits without any weights favour the hot reservoir. Therefore, the chosen weights ensure that our models properly fit all reservoirs, including the cold one that is best probed by the longer wavelengths that are more noisy given the lower sensitivity of JWST/MIRI’s Channel 4. While the choice for the weights is arbitrary, using lower weights will still favour the hotter reservoirs and result in a worse fit for the colder reservoir, whereas higher weights may start overrepresenting the cold reservoir and result in worse fits for the hot one.

In the MCMC, we used 20 walkers per free parameter and allowed the MCMC to explore the prior space for 50 000 iterations. Table 2 lists the prior spaces for each profile. Common priors that are kept the same in each fit are only listed once, for example, the number of molecules, N_mol, in the jump abundances. For some sources, we ran a second fit with a more limited prior space. In these cases, a second, less favourable solution was often found and a number of walkers got stuck in this local minimum, not allowing the fit to converge. The updated prior space was subsequently chosen to exclude these local minima. All fits use selected isolated H₂O lines and ortho-para line pairs, which all have been taken from Banzatti et al. (2025) (see their Figure 3 and Tables 5 and 6). The applied fit regions, which were also shifted by the heliocentric velocity, are listed in Table C.1. We used fewer regions for CX Tau and DN Tau, as the CO₂ P- and R-branches strongly contribute to their observed emission.

3.2 Rotational H₂O spectrum: Multiple components

As a comparison, we also followed Temmink et al. (2024a) to investigate the pure rotational H₂O spectrum using two (in the case of CX Tau and DN Tau) or three slab model components. Through trial and error (i.e. by fitting three or two components), we found that only two components were needed to fully describe the rotational H₂O spectra of CX Tau and DN Tau. This is very likely due to the H₂O emission being much weaker compared to the other sources (see Figure 1). We did not use all three approaches presented in Temmink et al. (2024a), but only the simplest one, approach I, which assumes a simple radial gradient. The line optical depths are sufficiently high that the vertical gradient in approach III does not matter and the results resemble those of approach I.

To agree with the parametric analyses discussed in Section 3.1, we slightly adapted approach I: instead of having a circular area surrounded by two annuli, we use three annuli with the inner radius of the first set to the dust sublimation radius, R_in. The other two or three radii are kept free, together with their respective temperatures and column densities. Therefore, Equation (2) of Temmink et al. (2024a) becomes: $\begin{array}{c}{F}_{\text{total }}=F_1\pi \left[{\left(\frac{R_1}{1\hspace{0.17em}\text{au}}\right)}^2-{\left(\frac{R_{\text{in}}}{1\hspace{0.17em}\text{au}}\right)}^2\right]\cos \hspace{0.17em}(i)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+F_2\pi \left[{\left(\frac{R_2}{1\hspace{0.17em}\text{au}}\right)}^2-{\left(\frac{R_1}{1\hspace{0.17em}\text{au}}\right)}^2\right]\cos \hspace{0.17em}(i)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+F_3\pi \left[{\left(\frac{R_3}{1\hspace{0.17em}\text{au}}\right)}^2-{\left(\frac{R_2}{1\hspace{0.17em}\text{au}}\right)}^2\right]\cos \hspace{0.17em}(i).\end{array}$(8)

For the sources for which we only fit two components, the third term in the above equation can be ignored.

The fits were done using the same MCMC setup, with the same reduced-χ² formula taken as the likelihood (see Equation (7)), as in Section 3.1. The same fit regions are used for the most optimal comparison between the results.

3.3 Line flux ratios

Aside from analysing the spectrum through slab models, we followed Banzatti et al. (2025) and analysed their suggested line flux ratios: F_{1500 K/3600 K} and F_{3600 K/6000 K} (see also the rightmost panel in their Figure 10). In particular, these line ratios will provide information on, respectively, the cold and hot H₂O vapour reservoirs. The temperatures indicate the approximate upper level energy. The 1500 K flux is comprised of two isolated H₂O transitions, with upper level energies of, respectively, 1448 K (at λ=23.81676 µm) and 1615 K (at λ=23.89518 µm), while the 3600 K and 6000 K fluxes are taken to be the isolated lines at 17.50436 µm and 17.32395 µm. Given their upper level energies, these lines probe the strength of the different temperature reservoirs and the ratios provide insights into the potential enhancement of the cold reservoir. Following Banzatti et al. (2025), we take the sum from the line fluxes of the 1448 K and 1615 K transitions as the total value for F_{1500 K}. To determine the line fluxes, we followed the approach of Banzatti et al. (2012) and Temmink et al. (2024b), based on the techniques provided by Pascucci et al. (2008), Najita et al. (2010), and Pontoppidan et al. (2010): we obtained a Gaussian distribution of measured line flux, using 1000 iterations of adding normally distributed noise (following the S /N values in Table B.1), by fitting Gaussians (using the python-package LMFIT, Newville et al. 2024) to the lines. From these distributions, we took the median value as the line flux and the full FWHM as the uncertainty.

4 Results

In the following section, we highlight the results from our different fitting methods. In Section 4.1, we present the parametric results for the different fit profiles, while the results of the multiple-component analysis are shown in Section 4.2. The calculated line fluxes for the different line tracers are displayed in Section 4.3.

4.1 Parametric analysis

Figure 3 displays all fitted profiles for our sample of compact discs with a line width of 4.71 km s⁻¹. Those for the quadrature line widths are given in Figure D.1. Given the profiles, we find three groups of discs: BP Tau, CY Tau, DR Tau, and RNO 90 all behave very similarly, which is most notable in the exponentially tapered power laws (second column), the jump abundances with T_jump at higher temperatures (T_jump ≳ 650 K, fourth column), and the resulting upward parabolas (fifth column). The resulting profiles for FT Tau and XX Cha are also nearly identical, but differ from the other four sources. The jump abundance occurs at a low temperature of T_jump ≤ 250 K for these discs, while their power laws are very slowly declining and their downward parabolas are very shallow. Finally, aside from their very similar spectra (see Figure 1), the parametrisations of CX Tau and DN Tau also behave very similarly to each other, but are different from the other discs. Their jump abundance occurs at the intermediate temperatures T_jump ∼ 500–600 K, while their power laws are monotonically increasing with radius and their parabolas are turned downwards.

The ${\chi }_{\text{red}}^2$ values are listed in Table D.1, where we have highlighted the parametrisations with the lowest ${\chi }_{\text{red}}^2$-value in boldface. Figures D.2–D.9 compare these best-fitting models with the observations and the fitting parameters are also listed in Tables D.2 and D.3. These tables also include uncertainty values on the retrieved fitting parameters. Those uncertainties are taken as the 16^th and 84^th percentiles of the posterior distributions and, therefore, represent the 1σ upper and lower limits, respectively. From the ${\chi }_{\text{red}}^2$ values, it is clear that not one single parametrisation provides the best description of all observations. Instead, a combination of the different profiles is preferred. We also note that the ${\chi }_{\text{red}}^2$ values of the different parametrisations are overall very similar.

Finally, we briefly highlight the resulting temperature slopes (values for q, see Equations (1) and (2)) according to our retrieved profiles. From Tables D.2 and D.3 it is clear that our retrieved values of q are between 0.35 and 1.30. Overall, we find the lower values for BP Tau, DR Tau, FT Tau, RNO 90, and XX Cha, while the higher values are found for CX Tau, CY Tau, and DN Tau. The average of the best-fitting parametrisations for all discs suggests a power law slope of q ∼0.57 for the temperature, which is somewhat lower than the value found by Romero-Mirza et al. (2024) of q ∼ 0.69 and the value suggested by thermochemical modelling of q ∼ 0.7 in the surface layers of discs (Bosman et al. 2022, Vlasblom et al., in prep.). Our average value would be even lower if we had only accounted for the results from the quadrature line widths.

4.2 Multiple component analysis

Overall, we find rather good agreement between the best-fitting profiles and the multiple components (three or two), as can be seen in Figures 3 and D.1, where the discrete data points indicate the results of such fits. The resulting values are given in Table D.4, while the ${\chi }_{\text{red}}^2$ values are also listed in Table D.1. We note that the radial locations shown in Figures 3 and D.1 are not the values for R₁, R₂, and R₃, the outer edge of each emitting area, but these are the central values of each emitting area. The resulting spectra for the best-fitting models are also included in Figures D.2–D.9 (red profiles).

The results for the multiple-component analysis are rather similar between the different discs with either three or two components fitted: for the three-component fits, we find that the temperatures for the first component all fall in the 765–970 K range, while those for the second and third components fall in the respective ranges of 355–530 K and 200–300 K. We note that the temperatures of the second and third components for FT Tau and XX Cha (T₂ ∼ 385 and T₃ < 215) are lower than those for BP Tau, CY Tau, DR Tau, and RNO 90, further suggesting a difference in behaviour as seen for the parametric analysis in Section 4.1. The temperatures for CX Tau and DN Tau, where only two components were fitted, are also very similar. The first component has a temperature between ∼ 400–500 K for both sources, with slightly higher temperatures for DN Tau, while the second component has a temperature in the range of 190–240 K.

Similarly, the column densities follow the same structure, where generally the first component has the highest column density, followed by the second and the third (where applicable). We note, however, that the retrieved column densities cannot be assumed to lie on a simple power law: as can be seen in Figure 3 and D.1, the other parametrisations provide better agreement with the multiple-component fits. For example, the downward parabola structure of CX Tau and DN Tau cannot be captured by a power law through the individual components. Finally, the high values for the emitting areas of the third component very likely indicate optically thin emission (see also Temmink et al. 2024a).

4.3 Line flux ratios

The values for the derived line fluxes are given in Table F.1, while the ratios are displayed in Figure 4, which has been adapted from Banzatti et al. (2025) (see their Figure 10, right panel) to include our sample of compact discs (the coloured data points). The grey data points are those from Banzatti et al. (2025), using the values listed in their Table 3. We find that BP Tau, CY Tau, DR Tau, and RNO 90 show similar line flux ratios as the sample of Banzatti et al. (2025), while the line fluxes of the cold lines (F_1500K) are much stronger in FT Tau and XX Cha. CX Tau and DN Tau also appear to have higher cold line fluxes, but we note that their hot H₂O line fluxes (F_6000K) are significantly lower. The implications of the locations of the different discs in Figure 4 are further discussed in Section 5.1.

Fig. 3 Profiles and multiple components (right-most panel) fitted for our sample of compact discs using a line width of 4.71 km s−1. The shaded area is the 1σ-confidence interval, given the uncertainties listed in Table D.2. The red stars in the top-left corner indicate the best fitting profiles for each disc (see Table D.1). We note that the best-fitting profile for DN Tau is highlighted in Figure D.1.

Fig. 4 Ratios of the 3600 K/6000 K and 1500 K/3600 K line fluxes, used to investigate the respective strength of each H2O reservoir. The grey data points are adapted from Banzatti et al. (2025) (see their Table 3), where we have highlighted those analysed by Romero-Mirza et al. (2024) and potential type P sources (i.e. MY Lup, RY Lup and IRAS-04385). The coloured data points indicate the line ratios of our sample of compact discs. See Section 5.1 for a discussion on the different types.

5 Discussion

In the following sections, we discuss and interpret our results. From our analysis, it immediately becomes clear that these discs, even though they are all classified as compact, show different behaviours given their H₂O spectra and inferred H₂O distributions. Not all of these discs show a strong, enhanced cold H₂O reservoir, suggesting that additional factors are at play within these systems. Within our sample of compact discs, many of them may behave just like large and structured discs, while only a limited number actually show cold H₂O emission to be strongly enhanced, likely due to radial drift.

5.1 Profile types and H₂O reservoirs

Given the different observed behaviours (see Section 4.1) and the observed line fluxes, we propose to categorise planet-forming discs in at least three different types, based on their H₂O reservoirs: type N, type E, and type P (as already highlighted in Figures 1 and 4). Type N includes the ‘normal’ discs, which show all three H₂O reservoirs (hot with T > 800 K, intermediate with 400 < T < 800 K, and cold with T < 400 K). The cold reservoir is present following the temperature gradient, expected to be present in all discs, but not strongly enhanced suggesting that radial drift is not efficiently replenishing the inner disc with icy pebbles. These discs can be found in the lower-left corner of Figure 4, where both ratios (F_{1500 K}/F_{3600 K} and F_{3600 K}/F_{6000 K}) have low values. We expect many of the large and structured discs to be part of this category.

The type E discs include those with very strong cold H₂O gas reservoirs, which appear to be significantly enhanced, suggesting efficient radial drift. These discs can be found in the top left corner of Figure 4, indicating their large cold H₂O reservoirs by large values for F_1500K/F_{3600 K}, while their values for F_3600K/F_6000K are similar to the type N discs. Finally, type P indicates the H₂O-poor (or CO₂-rich; Pontoppidan et al. 2010; Vlasblom et al. 2025) discs, whose spectra are generally devoid of H₂O, but still show strong emission from the cold reservoir. Emission from the other two components is often found to be weak and the emission present can be fitted with either a hot or an intermediate temperature. The lack of H₂O emission from both components (i.e. strong emission from both hot and intermediate) could potentially be due to a small inner cavity depleting the reservoirs in these components (Grant et al. 2023; Vlasblom et al. 2024). The cold component could still be prevalent due to radial drift, especially since snowlines are pushed further out in discs with cavities, holding for both the inner and outer regions (Temmink et al. 2023; Vlasblom et al. 2024). These H₂O-poor discs can be found in the top right corner of Figure 4, with high values for both F_1500K/F_3600K and F_3600K/F_6000K. We note that the boundaries for the different types in Figure 4 are arbitrary. To empirically determine the boundaries between the different types, a larger number of discs needs to be consistently analysed. We leave this for future work (Temmink et al., in prep.).

Based on our analysis of the 8 compact discs in the MINDS sample, we suggest that BP Tau, CY Tau, DR Tau, and RNO 90 can all be considered as type N discs, while FT Tau and XX Cha are the only two compact discs with very strong cold H₂O reservoirs: we consider them to be type E discs. CX Tau and DN Tau are considered to be type P discs. Using this categorisation we suggest that the type N discs can best be described (that is, lowest ${\chi }_{\text{red}}^2$-values) by profiles with jump abundances at higher temperatures or upward parabola, which also resembles the power law profile with a strong exponential taper. The type E discs are the ones with a jump abundance occurring at low temperatures (T < 250 K), have slow monotonically decreasing power laws, and may be described by a shallow downward parabola. On the other hand, the type P discs likely have rising power law profiles or jump abundances occurring at intermediate temperatures. Furthermore, their parabolas have a strong minimum.

The distinction between the different groups also becomes visible when investigating the contribution (with respect to the maximum flux in the observed spectrum) of the 50 individual slab models within the different profiles. As we are fitting 50 different slab models, if they all contributed equally, they would each reach a maximum contribution of 2% at every wavelength. Figure 5 displays the flux contribution of the slab models for DR Tau (top panel: type N), FT Tau (middle panel: type E), and CX Tau (bottom panel: type P). The contributions for DR Tau show that the hot component is the weakest (≤ 1.5%), as expected from the upward parabola, while the intermediate and cold reservoirs are of similar strength (∼ 3–4%) given the increasing emitting areas with decreasing temperature. For FT Tau, the hot and intermediate reservoirs are of similar strength (∼ 2%), while the cold reservoir has the strongest contribution (≤ 6%). The contributions for CX Tau show similar strengths in the hot and cold H₂O reservoirs (≤ 8%), but a negligible intermediate component, as expected from the downward parabola. While slab models probing the hot and cold reservoirs may have similar maximum contributions, fewer hot H₂O slabs contribute this much to the spectrum and the spectrum is dominated by the cold reservoir overall. Additionally, we want to emphasise that the best-fitting profile dictates the outcome of Figure 5: while the contributions for CX Tau, because of the parabola, have a strong contribution of the hot component, those of DN Tau (both type P discs) will lack the hot component and have significant contributions from the intermediate component as the spectrum is best-fitted by a jump abundance at a temperature of T_jump ∼ 600 K. Finally, we note that the line-to-continuum ratio of CX Tau is much lower (factor of ≥5) than for the other two sources, highlighting that even the higher percental contributions may not be as noticeable as for the other sources.

Fig. 5 Contributions of each individual slab model with the colour representing the excitation temperature of the respective slab. The contributions are shown for the best fitting parametric models of DR Tau (type N), FT Tau (type E), and CX Tau (type P). The contributions are taken as the percentage with respect to the maximum line flux in the full model. We note that the maximum line flux of CX Tau is at least a factor of 5 lower than that of the other two types.

5.2 Role of sub-structures

With our sample of compact discs now categorised, we highlight that the analysis of the ALMA visibilities of some of our sources has suggested that these rather compact discs may harbour sub-structures. In particular, gaps and rings have been proposed for the outer regions (>15 au) of BP Tau (Jennings et al. 2022; Zhang et al. 2023; Gasman et al. 2025), DN Tau (Long et al. 2018; Zhang et al. 2023), DR Tau (Jennings et al. 2022; Zhang et al. 2023; Gasman et al. 2025), and FT Tau (Long et al. 2018; Jennings et al. 2022; Zhang et al. 2023). With sub-structures hypothesised in discs across all three types, it is clear that the role of sub-structures in setting the inner disc H₂O reservoirs is not yet well understood (see also Gasman et al. 2025). Furthermore, as modelling works have shown that the influence depends on the gap location, the time at which the structure formed, and how leaky the traps are, this implies that sub-structures may play a very diverse role in setting the inner disc reservoirs (Kalyaan et al. 2023; Sellek et al. 2025; Easterwood et al. 2024). Age may also play a role (Mah et al. 2023), but we note that the majority of our sample consists of sources from the Taurus star-forming region and it can, therefore, be assumed that these discs all have a similar age. Therefore, our sample may be rather homogeneous and a larger sample, consisting of sources from a variety of different star-forming regions and with different ages, is needed to not only investigate the role of the structure formation time, but also the drift timescale (Mah et al. 2023, 2024).

Our introduced types may, therefore, be comprised of a wide variety of discs. Higher-resolution ALMA observations (spatial resolution of <0.1^′′, but preferably as high as possible) and observations with infrared observatories are needed to confirm and further study these sub-structures and their potential roles in setting the inner disc abundances, as well as hunt for sub-structures at smaller scales.

5.3 Profile implications

The following section discusses the differences in more detail. In particular, we further highlight the behaviours of the profiles and what this implies. In Appendix G, we also discuss how well the different profiles fit each type.

For the type N discs, we find that all discs prefer a jump abundance at higher temperatures of generally T_jump ≳ 700 K or a downward parabola (maxima generally at temperatures above 600 K). This behaviour suggests a potential depletion of the hotter component with respect to the stronger intermediate and cold components. This ‘depletion’ could simply mean that the hot H₂O reservoir is not as abundant as the intermediate and cold ones and, therefore, does not follow a simple power law profile, but it may also have a more physical explanation. As suggested by theoretical models (Sellek et al. 2025; Houge et al. 2025), the inward drift of icy pebbles not only increases the molecular reservoirs when crossing the snowlines, but also increases the continuum optical depth (τ) and, therefore, the continuum τ=1 layer may be moved to higher layers and obscure the underlying H₂O reservoirs.

On the other hand, the optical depth of the H₂O lines themselves may play an important role, which could reach values of τ ∼ 3000 (Temmink et al. 2024a). Therefore, it is possible that only the very top region of the emitting layer can be probed, well above the dust continuum τ=1 layer, while the remainder of the reservoir remains invisible given the high optical depths. As the line optical depth is found to be higher for the hottest component compared to the cooler components (Temmink et al. 2024a), this may play an important role in setting the observable column density for H₂O in the innermost region.

Finally, the destruction of H₂O could also play a role in depleting the inner regions. One of the main destruction routes in these high atmospheric layers of discs is photodissociation. The photodissociation must, however, occur on timescales faster than the gas-phase formation, which is efficient at >300 K (Bosman & Bergin 2021), to effectively destroy H₂O. Additionally, H₂O can be destroyed in these layers through collisions with C⁺, Si⁺, and H⁺ (Kamp et al. 2013). There is no consensus on which explanation for the flattening of the profiles in the inner regions is correct, as it may be a combination of all three: dust obscuration, optical depth, and destruction. We expect that it is more than likely that different combinations are preferred by different discs and the importance of each explanation will also differ.

The type E discs, FT Tau and XX Cha, have most notably the jump abundance at low temperature (T_jump < 250 K) or a shallow downward parabola. As drift can be expected to play to a certain extent a role in all discs, both smooth and structured, where the leakiness of the sub-structures determines the inward flow of material (Mah et al. 2023, 2024; Gasman et al. 2025), the cold component should be present in most if not all discs. The strength of the cold component in FT Tau and XX Cha, and therefore the jump in abundance at T_jump ∼ 250 K and the slowly decreasing power law, can be explained by an extreme or very efficient inward flux of the icy pebbles, enhancing the abundance of this cold component near the H₂O snowline with respect to the other discs.

We note that not all discs in the type E category may undergo efficient radial drift. Instead, an accretion outburst that significantly heats up the disc and, therefore, enlarges the emitting area of the cold component may also result in discs falling into this category. In particular, it has been found that XX Cha is highly variable in its accretion (Claes et al. 2022). Even though its variability timescales have not yet been constrained, the observed variation with the X-shooter instrument is on a scale of almost 2 dex (i.e. a factor of 100). Therefore, the strong cold component of XX Cha may be due to its strong variable accretion instead of efficient drift. Furthermore, recent work attributed the strong cold H₂O reservoir observed in the disc of EX Lup to an accretion outburst and the subsequent sublimation of H₂O ice, given the snowline being pushed outwards (Smith et al. 2025). For FT Tau, the other disc identified as a type E disc in our sample, there are currently no claims of accretion outburst or a strongly variable accretion rate. The reported accretion rate for FT Tau gives values between log₁₀ $\left({\dot{M}}_{\text{acc}}\right)\sim -7.5\hspace{0.17em}\text{and}\hspace{0.17em}-8.9$(Garufi et al. 2014; Gangi et al. 2022) (with ${\dot{M}}_{\text{acc}}$ given in M_⊙ year⁻¹), but have not been monitored closely over time. To confirm the role of accretion outbursts and variability in setting the inner disc chemical reservoirs, the accretion rates need to be monitored over an extended period of time to also infer the variability timescales. Also, if present, a larger sample of these sources needs to be observed with JWST-MIRI/MRS and analysed in a consistent manner.

The preferred parametrisations of CX Tau and DN Tau, the type P discs, include a monotonically increasing power law, a jump abundance at the intermediate temperatures, and a downward parabola. These profiles all imply a depletion of the hot or intermediate H₂O component with respect to the colder reservoir, which is clearly the strongest. This may suggest a potential dust and gas cavity (see, for example, Salyk et al. 2025), where the cold H₂O is most prominent due to sublimation of H₂O ice near the cavity edge (see also Vlasblom et al. 2024).

Finally, we note that Banzatti et al. (2023a) described two compact discs (GK Tau and HP Tau) having an ‘excess’ in their cold H₂O reservoir with respect to a large and structured disc (CI Tau). Here, we find that similar line flux ratios, as for their two compact discs and for our type N discs, do not necessarily imply an enhancement of the cold H₂O reservoir, but rather a smooth distribution from hot and/or warm to cold. A true enhancement, as found for our type E discs, produces an even stronger excess in the cold H₂O lines with respect to their warmer or hotter counterparts. This may also imply that CI Tau is actually depleted in the cold H₂O and, therefore, falls to the lower left corner of Figure 4. Therefore, the type E compact discs of FT Tau and XX Cha most closely conform with the hypothesis posed by Banzatti et al. (2020) and the overall situation appears to be more complex.

Overall, we find that the jump abundance with the jump temperature (T_jump) kept free is amongst the favoured fit for all three types. Therefore, we suggest that this profile can be fitted to investigate the H₂O reservoirs in all discs. Furthermore, the location of the jump, as discussed for each type above, may suggest which type the disc belongs to.

5.4 Line width comparison

One aspect that has not yet been addressed relates to the choices made for the different line widths, namely, between a fixed value of 4.71 km s⁻¹ or the sum in quadrature of a fixed turbulent line width and the thermal line width, given the temperature of the slab model. Intuitively, the quadrature line profiles may be preferred, as the line width changes accordingly with the temperature. One downside is that the strength of turbulence is not known in the inner regions of these discs, therefore, the turbulent component of the line width (fixed to 1 km s⁻¹) may be stronger or weaker. Thus far, a few studies have measured the turbulence in the outer regions of large discs using ALMA, resulting in values of <0.2 km s⁻¹ (Paneque-Carreño et al. 2024) or up to ∼ 35% of the sound speed (<0.35c_s; Flaherty et al. 2015, 2017, 2020). We find that seven out of eight of our discs prefer the fixed value of 4.71 km s⁻¹. This may suggest that the inner regions of these discs are more turbulent and higher values need to be assumed for the turbulent component (ΔV_turb > 1.0 km s⁻¹).

The preferred larger line width may be related to the change in the optical depth of the lines, which decreases with larger line widths (τ ∼ N/ΔV, with ΔV the line width). As the quadrature line widths are generally smaller (see also Romero-Mirza et al. 2024), these models will have larger optical depths. Furthermore, it becomes apparent from the contribution plots (Figures G.1–G.2) that the quadrature line profiles have stronger fluxes in certain lines of the cold component (see, for example, transitions around 25.05 or 26.70 µm). Given the added weights to the lines fitted at the longest wavelengths (>20 µm), we expect the fitting method to find the best profile describing the observed spectrum, regardless of the lower sensitivity of the larger wavelengths observed with JWST-MIRI. Therefore, it is possible that the cold component of the majority of our sample is less optically thick compared to the other discs.

Finally, we notice that our profiles show no qualitative differences between the different line widths (see Figures 3 and D.1), except for the parabola of DN Tau, which changes between upward and downward. Therefore, we conclude that the choice of line width does not really matter for the profiles, but note that the optical depth of the reservoirs, as discussed above, does depend on the chosen line width.

5.5 Comparisons with previous analyses

In this section, we compare our analysis with the works of Romero-Mirza et al. (2024) and Gasman et al. (2025), who, respectively, fitted (tapered) power-laws to seven other discs and analysed the H₂O emission in ten structured discs. Additionally, we compare our results with the interpretation of the models by Houge et al. (2025).

Romero-Mirza et al. (2024) only used power law and tapered power law profiles (our profiles I and II) to fit the column density distributions of four compact discs (FZ Tau, GK Tau, GQ Tau, and HP Tau) and three larger discs (AS 209, CI Tau, and IQ Tau). They find that all these discs can be well described by a simple power law, while both CI Tau and FZ Tau prefer the exponentially tapered power law, potentially suggesting a ring-like distribution of the H₂O reservoir. All discs are also captured in the analysis of Banzatti et al. (2025) and are, therefore, displayed in Figure 4, in which we have highlighted their respective positions. All of them, except for maybe GQ Lup, fall within our selection of type N discs. GQ Lup may also be a type N disc, given that it is located towards the lower right of Figure 4 with respect to FT Tau and XX Cha and has, therefore, a stronger intermediate component. As stated before, the current definition of the types on Figure 4 is somewhat arbitrary and refining the boundaries for each disc type requires a larger sample of sources, which will be addressed in future work (Temmink et al. in prep.). We also cannot exclude the possibility that other behaviours and types may exist outside the parameter space covered by the eight discs in our sample. As all discs analysed by Romero-Mirza et al. (2024) may be type N discs, their profiles can indeed be described by power laws (similar to CY Tau), but may also be represented by jump abundances at higher temperatures or upward parabolas.

Gasman et al. (2025) compared a sample of 10 structured discs and put their H₂O emission in context with respect to DR Tau. They note that the H₂O reservoirs are complicated for these structured discs and that the H₂O reservoirs strongly depend on the age of the system, the mechanisms that open a gap in the outer disc, and the leakiness of such sub-structures (see also Mah et al. 2024). We note that BP Tau and DR Tau are the two discs overlapping between their sample and ours, highlighting that compact discs may be structured as well. Relative to DR Tau, there is a distinct group of large, structured discs (CI Tau, DL Tau, GW Lup, and V1094 Sco) that have a weaker cold H₂O reservoir, which may end up in the lower left corner of Figure 4, where the position of CI Tau is already high-lighted. Another group of structured discs (SY Cha and Sz 98) shows, on the other hand, a very strong cold component (stronger than DR Tau) and may, therefore, end up in the upper left corner of Figure 4 near FT Tau and XX Cha. This suggests that the H₂O reservoirs of planet-forming discs are much more complicated and one cannot simply separate them into different types based on their size and/or dust structures, as already discussed in Section 5.3.

Finally, we compare our best-fitting profiles to the predicted model column density profiles by Houge et al. (2025) (for a star with L_∗=4 L_⊙, see Figure 6). Their column density profiles are given for the observable H₂O reservoir above the τ=1 layer at 20 µm. Their predictions are given for a fixed ratio between the chemical and mixing timescales (t_chem/t_mix=0.01) or for a fixed chemical timescale (t_chem=10 years) with a radially varying mixing timescale. Furthermore, their models include three types of dust models: fragile (v_frag=1 m s⁻¹), composition-dependent (v_frag=1–10 m s⁻¹), and resistant (v_frag=10 m s⁻¹). They find that for sufficient fragile dust (i.e. in the fragile or composition-dependent cases) the column density is set by the vapour-to-dust mass ratio and the dust opacity, as the profiles for those models do not strongly depend on any disc properties and, therefore, do not evolve over time (see also Sellek et al. 2025). While a direct comparison with the models of Houge et al. (2025) cannot be made given the stellar parameters used in their models, we compare the overall shape of the models with our best-fitting column density profiles in Figure 6. In particular, we make the comparisons with the models with a fixed chemical timescale and radially varying vertical mixing for the different dust types at a time of 1 Myr, an appropriate assumed age for the majority of our sources (see, for example, Krolikowski et al. 2021; Luhman 2023).

Qualitatively, the models and observations agree in the sense that the column density profiles decrease towards larger radii. However, our profiles decrease or flatten off towards the inner regions (see, for example, DR Tau and RNO 90). This could be due to the fact that the expression used by Houge et al. (2025) to derive the abundance of water in the surface layers assumes that the chemical timescale (i.e. the destruction timescale) is always much shorter than the mixing timescale. Under their assumptions, the abundance should saturate inside 0.16 au at a level equal to that of the midplane, resulting in the profile flattening off. This level, namely the limit of N ∼ 10²⁰ cm⁻² explored by Sellek et al. (2025), would be equivalent to negligible destruction under the assumption that the H₂O and dust particles behave the same and on similar timescales (that is, as long as the dust is not depleted on faster timescales than the H₂O). Our best-fitting profiles, except for the outer regions of CX Tau, do not exceed a column density of N ∼ 10²⁰ cm⁻², which is in agreement with the limit proposed by Sellek et al. (2025).

Quantitatively, the retrieved column densities fall below all the model profiles for many of our sources, even for the most conservative model predictions with the fragile dust. On the one hand, while their expected column density profiles are based on the increased dust obscuration following radial drift, the importance of the optical depth of the H₂O lines themselves is not explored, which may be important (as discussed in Section 5.3). Alternatively, this may suggest that for our sample of discs, the ratio of the chemical and mixing timescales is smaller than their assumptions. This may result from shorter photodissociation timescales, which are typically expected to be <1 year in the upper most layers of the disc (Bosman & Bergin 2021, Vlasblom et al. in prep.) and depend source-by-source on the stellar properties (i.e. the UV irradiation field). This may suggest that a chemical timescale of t_chem=10 yr may be too long. Future 2D models should further investigate the physical and chemical processing of H₂O in the inner regions and explore how effective the combination of increased dust obscuration following radial drift and line optical depth is.

Fig. 6 Comparisons between our parametric models with the lowest ${\chi }_{\text{red}}^2$-values and the prediction from the radially varying vertical models from Houge et al. (2025) for all three tested dust types: fragile (dotted), composition-dependent (solid), and resistant (dashed). The profiles of our discs identified to be type N are shown in purple, while those of the type E and type P discs are displayed in, respectively, blue and green.

5.6 Other molecular species

While our analysis focuses solely on the H₂O emission, the spectra of these compact discs contain more molecular species. In this section, we only highlight the other molecular species detected and leave the analysis of this emission for future work. We start by noting that CX Tau (Vlasblom et al. 2025) and DR Tau (Temmink et al. 2024b,a) have dedicated papers analysing their molecular emission and we, therefore, refer the reader to these papers. For the other discs in our sample, the commonly observed species (OH, CO₂, HCN, and C₂H₂) are detected in nearly all discs. For some discs, we also report the (potential) detection of CH₄ in CY Tau and that of some of the isotopologues, such as ¹³CCH₂ in CY Tau and both ¹³CO₂ and CO¹⁸O in DN Tau, which is very similar to CX Tau (Vlasblom et al. 2025). The different species detected in each disc are summarised in Table 3, while Appendix E contains a more elaborate discussion of the detected species.

6 Conclusions and summary

In this work, we analyse the pure rotational H₂O JWST-MIRI/MRS spectra of 8 millimetre-compact dust discs. We expand upon existing techniques to provide detailed parametric profiles of the column densities beyond (exponentially tapered) power laws. This leads us to infer the best-fitting radial profiles and investigate whether these compact discs show signatures of enhanced reservoirs due to radial drift or whether the overall situation is more complex. Our main conclusions are summarised as follows:

• The pure rotational H₂O spectra of compact dust discs are very different from each other. Half of the discs show similar strengths for the different H₂O reservoirs, conforming with what is seen among large and structured discs, while others show clear signs of an enhanced cold reservoir, even beyond that found in previous comparisons of compact and large discs;

• Different combinations of parametrisations can be used to describe the observed reservoirs in different discs, leading to the conclusion that planet-forming discs can be grouped into at least three different types (N, E, and P), based on their mid-infrared H₂O spectra;

• The parametrisation that assumes a constant of number of molecules with radius (i.e. the column density is proportional to negative the square of the radius) and has a jump (either an enhancement or depletion) in the abundance at a free temperature is able to provide a good fit for the column density profiles for all discs. The temperature of the jump varies per disc but is found to be similar within each of the different types (N, E, or P);

• We find that the column density profiles and, therefore, the distribution of the H₂O reservoirs are generally best fitted with a fixed line width of ΔV=4.71 km s⁻¹. As the quadrature sum is intuitively and physically more correct, this may suggest that the turbulence is stronger in the inner disc and, therefore, a larger value for the turbulent component needs to be assumed (ΔV_turb >1 km s⁻¹);

• Overall, a good agreement is found between the column densities retrieved through a multiple-component analysis and the profiles obtained through the parametric models. The parametric models show that one cannot assume a simple power law through the column densities derived from a multiple-component fit, but it is recommended to fit one of the other parametric models used throughout this work.

This work has shown that not all millimetre-compact discs have signatures of strongly enhanced cold H₂O reservoirs following radial drift. Instead, half of our sample behaves similarly to larger, more structured discs. Therefore, we have introduced a new categorisation, given the behaviour of the H₂O reservoirs seen for our discs: type N, type E, and type P discs. Here, N stands for ‘normal’ discs, whose spectra have contributions from all the components (hot, intermediate, and cold). Type P discs are the H₂O-poor discs, whose spectra are dominated by other molecular species, yet have a strong contribution from the cold reservoir. Finally, the type E discs are the ones that show enhanced cold H₂O, which may be the result of the strong inward drift of icy pebbles or, perhaps, the product of an accretion out-burst that pushes out the H₂O and increases the emitting area of the cold H₂O reservoir. Only two of our eight compact discs fit into this final category, further highlighting that not all compact discs have enhanced cold H₂O reservoirs. These three types provide a new classification of discs, purely based on their rotational H₂O emission and may be used to further explore trends involving system properties. Higher resolution ALMA observations and/or observations with infrared interferometers are required to further study the occurrence of sub-structures in these discs, on small and large scales, and how they may affect the observable reservoirs.

Data availability

Figures D.2–D.9 can be found on Zenodo: https://zenodo.org/records/15479557.

Acknowledgements

The authors would like to thank the referee for many thoughtful, constructive comments that helped improve the manuscript. This work is based on observations made with the NASA/ESA/CSA James Webb Space Tele-scope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1282. The following National and International Funding Agencies funded and supported the MIRI development: NASA; ESA; Belgian Science Policy Office (BEL-SPO); Centre Nationale d’Etudes Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luft- und Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economía y Competividad; Netherlands Research School for Astronomy (NOVA); Netherlands Organisation for Scientific Research (NWO); Science and Technology Facilities Council; Swiss Space Office; Swedish National Space Agency; and UK Space Agency. The data described here may be obtained from 10.17909/t6gq-q023. M.T., A.D.S., E.F.v.D., and M.V. all acknowledge support from the ERC grant 101019751 MOLDISK. D.G. thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA). E.F.v.D. also acknowledges support the Danish National Research Foundation through the Center of Excellence “InterCat” (DNRF150). E.F.v.D., I.K., and A.M.A. acknowledge support from grant TOP-1 614.001.751 from the Dutch Research Council (NWO). T.H., K.S. and M.S. acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. A.C.G. acknowledges support from PRIN-MUR 2022 20228JPA3A “The path to star and planet formation in the JWST era (PATH)” funded by NextGeneration EU and by INAF-GoG 2022 “NIR-dark Accretion Outbursts in Massive Young stellar objects (NAOMY)” and Large Grant INAF 2022 “YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems (YODA)”. G.P. gratefully acknowledges support from the Carlsberg Foundation, grant CF23-0481 and from the Max Planck Society. B.T. is a Laureate of the Paris Region fellowship program, which is supported by the Ile-de-France Region and has received funding under the Horizon 2020 innovation framework program and Marie Sklodowska-Curie grant agreement No. 945298. This work also has made use of the following software packags that have not been mentioned in the main text: NumPy, SciPy, Astropy, Matplotlib, pandas, IPython, Jupyter (Harris et al. 2020; Virtanen et al. 2020; Astropy Collaboration 2013, 2018, 2022; Hunter 2007; pandas development team 2020; Pérez & Granger 2007; Kluyver et al. 2016).

Appendix A Velocity shifted slab models

Fig. A.1 Comparison, using DR Tau as a test case, between an unshifted slab model (red) and one shifted by the heliocentric velocity.

Appendix B S/N on the continuum and uncertainties for the observations

Appendix C Slab model fit regions

Appendix D Fit results

Fig. D.1 Similar to Figure 3, but for the quadrature line widths. The shaded area is the 1σ-confidence interval, given the uncertainties listed in Tables D.3.

Appendix E Other molecular species

In this following section, we discuss the observed molecular species, aside from H₂O, in each disc (Section E.1) and we discuss the tentative detection of CH₄ in CY Tau (Section E.2). These features are also highlighted in Figures D.2-D.9.

E.1 Molecules per source

BP Tau Aside from the pure rotational H₂O lines, the spectrum of BP Tau also contains strong emission of the ro-vibrational transitions and CO emission at the shortest wavelengths (4.9-5.3 µm). Other molecular species, such as CO₂ and HCN, are only weakly detected, while OH transitions are strongly detected above >15 µm. Both the H₂ 0-0 S(1) (at 17.30484 µm) and S(2) (at 12.27861 µm) are detected. The detected molecular species at >10 µm are also highlighted in Figure D.2.

CY Tau The strongest molecular feature in the spectrum of CY Tau belongs to C₂H₂, also visible in Figure 1. Its isotopologue, ¹³CCH₂ is also well detected. Other molecular species include CO, the ro-vibrational lines of H₂O, OH, HCN, and CO₂. Both the H₂ S(1) and S(2) transitions are also detected. Additionally, we present a potential detection of CH₄ in this disc. The tentative detection is further highlighted in Section E.2. The molecular species are highlighted in Figure D.4.

DN Tau DN Tau has previously been classified, just as CX Tau, as a CO₂-rich source (Pontoppidan et al. 2010). Aside from the strong emission of CO₂, we report the detection of its isotopologue ¹³CO₂ and the potential detection of CO¹⁸O. Other molecular species, such as HCN and C₂H₂ are not confidently detected, while OH is. At the shortest wavelengths, CO is only tentatively detected, while emission from the ro-vibrational H₂O transitions is clearly detected. As for the other sources, both the H₂ S(1) and S(2) transitions are also detected. Figure D.5 high-lights the detected molecular species.

FT Tau Aside from the strong pure rotational H₂O emission, the spectrum of FT Tau also contains emission from the ro-vibrational lines and from CO. At the longer wavelengths, we clearly detect OH, CO₂, HCN, and C₂H₂. The H₂ S(1) and S(2) transitions are also detected. The emission of these species is also highlighted in Figure D.7.

RNO 90 Besides the rotational H₂O transitions, the spectrum of RNO 90 also contains emission from CO, the rovibrational H₂O transitions, OH, CO₂, HCN, and C₂H₂. Both the S(1) and S(2) transitions of H₂ are also detected. The molecular emission is highlighted in Figure D.8.

XX Cha Although the spectrum of XX Cha is very similar to that of FT Tau, the emission from both CO₂ and C₂H₂ is stronger. Emission from CO, the ro-vibrational H₂O lines, OH, and HCN are also clearly detected. As for all the other sources, both the H₂ S(1) and S(2) transitions are also detected. Figure D.9 highlights the molecular emission of the other species in the spectrum of XX Cha.

E.2 Tentative detection of CH₄ in CY Tau

We report a tentative detection of CH₄ in the spectrum of CY Tau, whose Q-branch is located at ∼ 7.66 µm. In Figure E.1, we display the 7.60-8.00 µm wavelength region, containing molecular emission from ro-vibrational H₂O transitions and potential weak features from C₂H₂ and the potential detection of CH₄. For the slab models (H₂O, C₂H₂, and CH₄) we have used excitation temperatures of 975 K, 500 K, and 500 K, column densities (log₁₀(N) with N in cm⁻²) of 18.4, 17.0, and 16.8, and emitting radii of 0.03 au, 0.10 au, and 0.10 au, respectively. The parameters for H₂O and C₂H₂ were obtained using a ${\chi }_{\text{red}}^2$ approach similar to that described in Section 3.1 but using grids instead of an MCMC exploration (see also, for example, Grant et al. 2023). The model for CH₄ has been obtained by visually fitting the residuals. These slab models thus suggest a column density ratio of $N_{{\text{C}}_2{\text{H}}_2}/N_{{\text{CH}}_4}\gtrsim 1.5$, under the assumption that the species are co-existing. The potential presence of CH₄, in addition to the detection of ¹³CCH₂, suggests that CY Tau may be the most carbon-rich disc in our sample. This disc may be similar to the much larger disc of DoAr 33, where Colmenares et al. (2024) reported, alongside other hydrocarbons, a tentative detection of CH₄. They compared their spectra to thermochemical models and argued for an overall carbon-to-oxygen ratio (C/O) larger than unity.

Fig. E.1 Zoom-in on the 7.60-8.00 µm region of CY Tau, showing the potential detection of CH4 (green). Also shown are approximate contributions from H2O (blue) and C2H2 (yellow) slab models. The full model is shown in red with a small offset applied to the flux.

Appendix F Integrated line fluxes

Appendix G Profile comparison

In Figures G.1, G.2, and G.3, we show the flux contributions of the 50 slab models for the profiles fitted for DR Tau, CX Tau, and FT Tau, respectively. These discs are each part of one of the different types and comparing their contribution plots provides information on how well the different profiles fit for their types.

DR Tau (type N) The type N discs (BP Tau, CY Tau, DR Tau, and RNO 90) are overall best described by the exponentially tapered power laws, the jump abundances with the jump occurring at high temperatures, and upward parabola (see also Table D.1). As can be seen in Figure G.1, those profiles all yield very similar contributions. The main differences between those profiles and the ones with the simple power law or the jump abundance with the jump fixed at T_jump=400 K is the flux at the shortest wavelengths (<15.0 µm), best probed by the hottest component, where the power law and the fixed temperature jump abundance have stronger contributions of the hotter components (most notable in the left panels of Figure G.1. These stronger contributions results in overfitting the observed flux in those inner regions, whereas the preferred profiles have a decrease in column density in the innermost and, therefore, do not overfit these inner regions.

FT Tau (Type E) FT Tau and XX Cha, the Type E discs, are best described by a jump abundance at low temperatures or parabola with a small upturn (see also Table D.1). As can be seen in Figure G.2 for FT Tau, the jump abundance with a jump at T_jump ≲250 K is really able to capture the strength of the cold component. The power law is not able to capture this cold component as strongly, but the contributions are clearly still there. The parabola, while being able to yield a similar contribution to the cold component as the power law, has a much stronger contribution at the shortest wavelengths (<15.0 µm) following its upward turn. This stronger contribution actually overfits the observed flux at these wavelengths, disfavouring the parabola as the best fit.

CX Tau (type P) The type P discs (CX Tau and DN Tau) mainly prefer the downward parabola and jump abundance profiles (see also Table D.1). As can be seen in Figure G.3, those profiles have the strongest contributions from either the hot or intermediate and cold components, while the (tapered) power law profiles are not able to capture the strength of all three components (hot, intermediate, and cold). The preference for these profiles shows that the P type discs do have contributions from at least two of the components.

Fig. G.1 Contributions from all parametric fits for DR Tau, representing the type N discs.

Fig. G.2 Similar as Figure G.1, but for FT Tau, representing the type E discs.

Fig. G.3 Similar as Figure G.1, but for CX Tau, representing the type P discs.

References

All Tables

All Figures

	Fig. 1 Normalised-to-peak-flux spectra of our sample of millimetre-compact discs. The red line indicates the estimated continuum. The indicated types and their meaning are discussed in Section 5.1. In blue, we highlight the 23.72–24.03 µm wavelength region, where two transitions are located that are most important for obtaining information about the cold H2O reservoir.
In the text

	Fig. 2 Schematic highlighting the different parametric models used for the column density profiles. Both the horizontal (radius) and vertical (column density) axes are taken in log10-space.
In the text

Fig. 3 Profiles and multiple components (right-most panel) fitted for our sample of compact discs using a line width of 4.71 km s−1. The shaded area is the 1σ-confidence interval, given the uncertainties listed in Table D.2. The red stars in the top-left corner indicate the best fitting profiles for each disc (see Table D.1). We note that the best-fitting profile for DN Tau is highlighted in Figure D.1.

In the text

Fig. 4 Ratios of the 3600 K/6000 K and 1500 K/3600 K line fluxes, used to investigate the respective strength of each H2O reservoir. The grey data points are adapted from Banzatti et al. (2025) (see their Table 3), where we have highlighted those analysed by Romero-Mirza et al. (2024) and potential type P sources (i.e. MY Lup, RY Lup and IRAS-04385). The coloured data points indicate the line ratios of our sample of compact discs. See Section 5.1 for a discussion on the different types.

In the text

Fig. 5 Contributions of each individual slab model with the colour representing the excitation temperature of the respective slab. The contributions are shown for the best fitting parametric models of DR Tau (type N), FT Tau (type E), and CX Tau (type P). The contributions are taken as the percentage with respect to the maximum line flux in the full model. We note that the maximum line flux of CX Tau is at least a factor of 5 lower than that of the other two types.

In the text

Fig. 6 Comparisons between our parametric models with the lowest ${\chi }_{\text{red}}^2$-values and the prediction from the radially varying vertical models from Houge et al. (2025) for all three tested dust types: fragile (dotted), composition-dependent (solid), and resistant (dashed). The profiles of our discs identified to be type N are shown in purple, while those of the type E and type P discs are displayed in, respectively, blue and green.

In the text

	Fig. A.1 Comparison, using DR Tau as a test case, between an unshifted slab model (red) and one shifted by the heliocentric velocity.
In the text

	Fig. D.1 Similar to Figure 3, but for the quadrature line widths. The shaded area is the 1σ-confidence interval, given the uncertainties listed in Tables D.3.
In the text

	Fig. E.1 Zoom-in on the 7.60-8.00 µm region of CY Tau, showing the potential detection of CH4 (green). Also shown are approximate contributions from H2O (blue) and C2H2 (yellow) slab models. The full model is shown in red with a small offset applied to the flux.
In the text

	Fig. G.1 Contributions from all parametric fits for DR Tau, representing the type N discs.
In the text

	Fig. G.2 Similar as Figure G.1, but for FT Tau, representing the type E discs.
In the text

	Fig. G.3 Similar as Figure G.1, but for CX Tau, representing the type P discs.
In the text