© The Authors 2025

1 Introduction

Understanding the chemistry of the inner planet-forming region of protoplanetary disks is crucial for understanding the composition and atmospheric properties of exoplanets. Most rocky and super-Earth planets are believed to form within the inner few au of disks and thereby inherit their bulk composition directly from this region (e.g., Dawson & Johnson 2018; Öberg & Bergin 2021). The inner disk regions can be probed with infrared (IR) facilities and are now characterized in greater detail by the James Webb Space Telescope (JWST) (e.g., Kamp et al. 2023; Henning et al. 2024; van Dishoeck et al. 2023; Banzatti et al. 2023; Romero-Mirza et al. 2024; Gasman et al. 2025; Arulanantham et al. 2025). In particular, the characterization of H₂O emission with JWST is of great interest due to its potential role in creating habitable conditions (see, e.g., Krijt et al. 2023).

The rotational H₂O spectrum – as seen with, for instance, JWST/MIRI – covers transitions spanning a wide range of upper level energies E_up (~1000–10 000 K) and Einstein-A coefficients A_ij (~0.01–100 s⁻¹). These lines hence probe a wide range of excitation conditions (see, e.g., Gasman et al. 2023; Banzatti et al. 2023, 2025; Temmink et al. 2024). Using the relative excitation of higher and lower E_up lines, the H₂O spectrum can be used as a proxy to retrieve the temperature and H₂O column density structure of the IR emitting region in the inner disk (Romero-Mirza et al. 2024; Temmink et al. 2025). It thus contains very valuable information, both for informing disk modeling efforts and for understanding the conditions under which nascent planets may be forming in the disk.

H₂O is produced in the inner disk via high-temperature gasphase formation at rapid timescales via the OH + H₂ ⟶ H₂O + H reaction, which is efficient at temperatures ≳300 K (Charnley 1997; van Dishoeck et al. 2013; Walsh et al. 2015). In the warm, upper layers of the disk, this reaction is balanced by rapid photodissociation due to the lower dust opacity. In the outer disk (greater than a few astronomical units), the temperature becomes low enough (≲150 K) for the H₂O to freeze out onto the dust grains. Above this ice reservoir, only a cold photodesorption layer is present, as the temperature is too low in this region for either gas-phase formation or thermal desorption (Woitke et al. 2009; Hogerheijde et al. 2011; van Dishoeck et al. 2021). Aside from these chemical processes, the H₂O vapor in the inner disk is also affected by dynamical processes, such as the inward drift of dust. When the icy grains in the outer disk move inward, their ice mantles can sublimate once the grain crosses a species snow line. This can enhance the amount of H₂O vapor in the inner disk and has been explored using data from the Spitzer Space Telescope by Banzatti et al. (2020), who find a correlation between the H₂O flux and dust disk size. Disks with a smaller size as measured by the millimeter dust continuum, thought to be caused by efficient radial drift, are found to have stronger H₂O fluxes, suggesting a link between dust drift and H₂O vapor in the inner disk. This correlation has also been theoretically motivated by Kalyaan et al. (2021). The effects of dust drift and H₂O delivery to the inner disk have been further explored in 1D by the modeling work of, for example, Mah et al. (2023), Sellek et al. (2025), and Houge et al. (2025), and the effects of dust evolution on mid-IR lines has also been explored in 2D by, for example, Greenwood et al. (2019).

More recently, this link between dust dynamics and H₂O emission has been expanded upon with the help of JWST’s improved spectral resolution and sensitivity compared with Spitzer. Banzatti et al. (2023) and Romero-Mirza et al. (2024) demonstrate that compact dust disks exhibit excess H₂O flux in the cold, low-E_up lines with respect to more extended dust disks and propose that this is a signature of strong drift. The low-E_up lines trace temperatures of ≲200–250 K (see also Banzatti et al. 2025), which likely trace close to the H₂O snow line. A large volume of ice crossing the snow line may thus lead to an increased amount of cold, ~200 K H₂O vapor, which translates to an increase in cold H₂O flux.

Most work so far has used local thermodynamic equilibrium (LTE) slab models to characterize the H₂O emission in JWST/MIRI spectra. These models parameterize the emission using only a temperature (T), column density (N), and emitting area (A). They can be used in various ways such as single-component fits (e.g., Grant et al. 2023; Tabone et al. 2023) and two- or three-component fits (e.g., Temmink et al. 2024; Pontoppidan et al. 2024). As these models are used to inform us of the temperature and column densities within the planet-forming zone, understanding the accuracy of these retrievals is key to our understanding of disk chemistry.

2D thermochemical models provide an ideal case to put these retrieval methods to the test. The retrieved temperature and column density provide at most only a 1D view of the disk and thus do not directly probe the underlying 2D distribution of H₂O vapor in the disk. Thermochemical models do contain the full 2D view of the disk, thereby providing insights that slab models alone cannot give. They can inform us how retrievals may provide indirect information about the underlying 2D abundance structure, and they can inform us of the accuracy of the retrievals by directly comparing the retrieved temperatures (the excitation temperature, T_ex) and column densities to the model output (e.g., the kinetic temperature, T_kin). As a further test, we also investigated the retrievals on CO₂ emission (see Appendix A.1).

The goal of this paper is two-fold: we aim to provide a better understanding of the accuracy of retrieval techniques commonly used on observations and also investigate the cold H₂O lines that have been proposed as signatures of radial drift, in order to understand the insights that thermochemical models can provide on the origin of this emission. The structure of this work is as follows: in Sect. 2, we describe our modeling approach using the Dust And LInes (DALI) thermochemical code, as well as the retrieval techniques we used to interpret the synthetic DALI spectra. In Sect. 3, we describe the results of our modeling, retrievals, and our findings on the cold H₂O emission. We discuss the origin of this cold H₂O and its potential relation to radial drift further in Sect. 4 and also provide predictions for observing cold H₂O lines with future far-IR missions. We summarize our main conclusions in Sect. 5.

2 Methods

2.1 Thermochemical models

This work makes use of the 2D thermochemical code DALI (see Bruderer et al. 2012; Bruderer 2013). This code consists of three main steps. First, the code calculates the local radiation field and dust temperature at all locations in the disk given an input gas and dust density structure. Second, the gas temperature, chemical abundances of all species in the chemical network, and non-LTE excitation for several specified atoms and molecules (and thus their cooling ability) are calculated self-consistently. Third, a ray tracing tool is available to obtain line fluxes, spectra and spectral image cubes. We present two sets of models: in the first set of models, the chemical abundances of all species are set using the chemical network and thermal balance; in the second set, the 2D abundances are parameterized, as this allows us to investigate more directly how differences in the abundance structure manifest themselves in the spectrum.

The input parameters of all models are summarized in Table A.2. They are based on those presented in Vlasblom et al. (2024), so we refer the reader to that work and references therein for the full details of the model setup. These models are based on the model by Zhang et al. (2021) and Bosman et al. (2022a,b) for the disk around K5 star AS 209 and therefore use the input stellar spectrum from those works, which has an effective temperature of 4300 K and a bolometric luminosity of 1.4 L_⊙. The models contain two dust populations, separated by their grain size. The small dust population has sizes between 5 nm and 1 μm, and is assumed to be well-coupled to the gas, following the same radial and vertical distribution (Miotello et al. 2016). The large dust population has sizes between 5 nm and 1 mm and is assumed to be more settled to the midplane than the small dust, and thus has a reduced scale height. Both populations follow an MRN distribution (Mathis et al. 1977). The large grains have a mass fraction f_ℓ, which sets the gas-to-dust ratio in the upper layers of the disk due to the settling of the large grains. The models presented in the main body of this work have f_ℓ = 0.9, which produces a gas-to-dust ratio of ~10³ in the IR emitting region of the disk. This is in line with the predicted elevated gas-to-dust ratio in these layers from Meijerink et al. (2009) and Woitke et al. (2018). In Appendix A.3, we discuss the effects of even higher gas-to-dust ratios using models with f_ℓ = 0.99 and 0.999. For our first set of models where we run the chemistry and thermal balance, we use the same chemical network as Bosman et al. (2022a,b). This is a network that builds on the standard DALI chemical network – based on UMIST06 (Woodall et al. 2007) – by including the effects of H₂O UV-shielding (Bethell & Bergin 2009) and more efficient H₂ formation at high temperatures. We note that the use of a more recent UMIST version, such as UMIST12 or UMIST22, will not affect the results of this work. Isotope chemistry (as implemented in, e.g., Miotello et al. 2014; Visser et al. 2018) is not included in this network.

Aside from the set of models containing the full chemistry, we also present a second set of models where the abundance structure is parameterized using a jump abundance profile, following the methods in Bruderer et al. (2015) and Bosman et al. (2017). In these models, the second step of the code (the chemistry and gas thermal balance) is thus skipped, meaning that the gas temperature cannot be calculated self-consistently. Prior studies have typically assumed it to be equal to the dust temperature. However, in the IR emitting layer, the gas and dust are thermally decoupled, making this assumption invalid (see, e.g., Bosman et al. 2022a; Vlasblom et al. 2024). Thus, we used the self-consistently calculated temperature profile from our full chemistry models as the assumed temperature structure for our parameterized models. This assumption gives a much more realistic gas temperature than setting T_gas = T_dust. However, we do note that our assumed gas temperature for the parameterized models is not entirely self-consistent, as any change in abundance structure could slightly alter the gas temperature as well due to H₂O being both an effective coolant as well as a source of heating through its photodissociation in the surface layers (see, e.g., Woitke et al. 2018, 2024).

The parameterized abundance structures are inspired by the chemical network. We parameterized the H₂O abundance structure using two reservoirs. The first is an inner reservoir with a high fractional abundance of n(H₂O)/n_H = 10⁻⁴ (where n_H = n(H) + 2n(H₂)) defined as the region where A_V > 1.5 mag and T_dust > 150 K. In the rest of the disk, the fractional H₂O abundance is set to 10⁻¹⁰. In Appendix A.1, we show how we parameterized the CO₂ abundance structure. The abundance structure of H₂O for both the parameterized and full chemistry models are presented in the left column of Fig. 1.

We made use of DALI’s “fast ray-tracer” (see Bosman et al. 2017, Appendix B) to generate the IR spectra of H₂O (as well as ¹²CO₂ and ¹³CO₂) between 10 and 28 μm, assuming a distance of 121 pc and a face-on orientation. Since isotope chemistry is not included in our chemical network, the ¹³CO₂ spectrum was derived by scaling the ¹²CO₂ abundance down by a factor of 77, the interstellar medium (ISM) ¹²C/¹³C ratio (Wilson & Rood 1994). The excitation of the relevant molecules was calculated in non-LTE for both sets of models. The molecular data file for H₂O was retrieved from the Leiden Atomic and Molecular Database (LAMDA), including levels with energies up to 7200 K (Tennyson et al. 2001). The line transitions were obtained from the BT2 list (Barber et al. 2006) and the collisional rate coefficients come from Faure & Josselin (2008). For CO₂, we used the molecular data compiled by Bosman et al. (2017) – who retrieve the energy levels, line positions, and line strengths from the HITRAN database (Rothman et al. 2013) – and used collisional rate coefficients based on Allen et al. (1980), Nevdakh et al. (2003), and Jacobs et al. (1975). We compared the non-LTE results with those of an LTE excitation calculation in Appendix A.2. We convolved the generated spectra to the approximate resolving power of MIRI of R = 3000, and we did not add synthetic noise.

Fig. 1 Abundance maps of H2O (left and middle panels) and temperature map (right panel) for the fiducial models with fℓ = 0.9. The left panel depicts the model using the chemical network, and the middle depicts the model with parameterized abundances. In all panels, the dust τ = 1 surface at 15 μm and the AV = 1.5 surface are indicated with a dashed blue line and a dotted black line, respectively. The red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371 K), 113,9−100,10 (Eup = 2438 K), and 83,6−70,7 (Eup = 1447 K) lines, respectively, representing hot, warm, and cold H2O. The warm H2O emitting region for both models is also depicted in the right panel in black. The light blue line in the left and middle panels represents the H2O snow surface.

2.2 Slab retrievals

Having generated our synthetic spectra, we performed several retrieval analyses using LTE slab models to better understand what those models trace in terms of regions and conditions. These slab models contain three free parameters: the gas temperature (T), the column density (N), and the emitting area (A) (see, e.g., Kamp et al. 2023). If one assumes $A=\pi R_{\mathrm{e}\mathrm{q}}^2$, the emitting area can be converted to an equivalent slab radius (R_eq). This equivalent radius can reflect the true emitting radius of the emission, but it can also originate from an annulus further out in the disk instead. We assumed Gaussian line profiles with a full width at half maximum of ΔV = 4.7 km s⁻¹, following Salyk et al. (2011). In the main body of this paper, we focus on our analysis of the H₂O emission. We exclude the ro-vibrational band at ~6–8 μm from our analysis as this band is known to be much more affected by non-LTE effects (see, e.g., Meijerink et al. 2009; Pontoppidan et al. 2024; Banzatti et al. 2025). In Appendix A.1, we present our analysis of the CO₂ emission.

The H₂O rotational spectrum contains a plethora of lines with a wide range in upper level energy, E_up, thereby tracing a wide range of temperatures and densities (see, e.g., Banzatti et al. 2023, 2025). Thus, fitting the entire rotational spectrum from 10 to 28 μm with a single temperature is unlikely to provide the best insight. Hence, we used two different approaches.

First, we performed single-temperature fits in three different wavelength ranges (10–14, 13.5–17.5, and 21–24 μm), following studies such as Gasman et al. (2023) and Temmink et al. (2024). As seen in, for instance, Banzatti et al. (2023), the H₂O lines at the shortest wavelengths (~10 μm) generally have the highest E_up and should thus trace the highest temperatures in the innermost regions of the disk, with subsequent lines at longer wavelengths having lower E_up and thus tracing further out. The best-fit model was obtained from a grid of models using a χ² method. In this grid, T was varied linearly in 40 intervals between 100 and 1000 K (ΔT = 22.5 K) and N was varied in log-space in 40 intervals between 10¹⁴ and 10²³ cm⁻² (Δ log N = 0.225). The spacing of this grid was kept consistent for all fits. For each point in the grid, the best-fit R_eq was determined by minimizing the reduced χ². We did not account for the effects of mutual shielding of adjacent lines (as in, e.g., Tabone et al. 2023), as this was also not accounted for when generating our synthetic DALI spectra. Banzatti et al. (2025) demonstrate that, for H₂O, mutual shielding is important for individual orthopara line pairs (see their Fig. 5) but does not significantly affect the spectrum as a whole. Temmink et al. (2024) also demonstrate that retrievals are only negligibly affected.

Second, we fit the entire rotational spectrum using three temperature components following scenario I presented in Temmink et al. (2024). This method combines the flux from three different slab models with decreasing temperature and increasing emitting area, where the best fit is found using the Monte Carlo Markov chain (MCMC) implementation emcee (Foreman-Mackey et al. 2013). The total H₂O flux is parameterized as $$\begin{array}{rl}{F}_{\text{total\hspace{0.17em}}}=F_1\pi {\left(\frac{R_1}{1\mathrm{a}\mathrm{u}}\right)}^2& +F_2\pi [{\left(\frac{R_2}{1\mathrm{a}\mathrm{u}}\right)}^2-{\left(\frac{R_1}{1\mathrm{a}\mathrm{u}}\right)}^2],\\ & +F_3\pi [{\left(\frac{R_3}{1\mathrm{a}\mathrm{u}}\right)}^2-{\left(\frac{R_2}{1\mathrm{a}\mathrm{u}}\right)}^2],\end{array}$$(1)

where F_i and R_i represent the flux and emitting radii of the three slab model components. We did not take into account the effects of spatial shielding of the colder components by the hotter ones (Scenario III in Temmink et al. 2024), as this was found to be negligible by Temmink et al. (2024) due to the high optical depth of the hot components.

Fig. 2 Synthetic H2O (blue), 12CO2 (green), and 13CO2 (purple) spectra. The top row depicts the model using the chemical network, and the bottom row depicts the model with parameterized abundances. In the right panels, blue triangles indicate H2O lines with Eup < 2500 K.

3 Results

3.1 Abundance maps and synthetic spectra

We present the H₂O abundance maps for our full chemistry (left panel) and parameterized abundance models (middle panel), as well as the gas temperature profile (right panel) in Fig. 1. The red, pink, and blue contours in the first two panels indicate the line emitting region from which 70% of the emission originates. These contribution functions (CBFs; see Bruderer 2013 for details on its calculation) are shown for three H₂O lines with different E_up, namely the 17_7,10−16_4,13 (6371 K; red), 11_3,9−10_0,10 (2438 K; pink), and 8_3,6−7_0,7 (1447 K; blue) lines (see Table A.1 for an overview of the properties of all analyzed lines in this work). These lines trace the warm surface layer of the disk, above the τ = 1 surface of the dust continuum at 15 μm. As expected, the lines with higher E_up trace hotter regions closer to the star, whereas the lower-E_up lines trace slightly further out in the disk, but still within 1 au.

The synthetic spectra for both sets of models are presented in Fig. 2. In the left panel, the H₂O, ¹²CO₂, and ¹³CO₂ spectra are shown between 13.75 and 16.5 μm. When comparing the full chemistry and parameterized abundance models, the flux ratios between the three species are similar. The CO₂ emission is stronger than that of H₂O, which can likely be attributed to the gas-to-dust ratio in the surface layers of these models, which is set to ~10³. Vlasblom et al. (2024) demonstrate that higher gas-to-dust ratios (~10⁴−10⁵) increase the temperature in the surface layers (see also Woitke et al. 2018), as well as the H₂O/CO₂ flux ratio, likely due to the higher temperatures promoting H₂O production over CO₂.

Interestingly, when considering the H₂O emission between 21 and 24 μm (right panels of Fig. 2), the full chemistry and parameterized abundance models are less similar. In the parameterized model, the cold H₂O emission (lines with E_up < 2500 K are marked by blue triangles) is strong compared to the surrounding H₂O lines with higher E_up. Moreover, the H₂O lines at ~15 μm (which are generally dominated by higher E_up than those at ~23 μm) are much weaker than the lines at longer wavelengths, which indicates that the parameterized model is dominated by strong cold H₂O. In the full chemistry model, the line fluxes at ~15 and ~23 μm are much more similar and are indicative of warmer emission. This must be related to the difference in the abundance structure between the two models, as all other factors such as the temperature structure are controlled for. We discuss this in detail in Sect. 3.3.

3.2 Slab retrievals on synthetic spectra

3.2.1 Temperature

To compare the results from the slab retrievals to our thermochemical models, we first generated a 1D temperature profile from our 2D temperature structure using the CBFs (see the contours plotted on the abundance maps in Fig. 1). We extracted the 1D temperature profile by calculating the CBF-weighted average temperature within the 70% CBF at each radius, following Kaeufer et al. (2024). This was done for all three aforementioned lines. As these lines generally trace similar depths and their radial ranges partially overlap, the three derived temperature profiles closely agree with one another. We compare these temperature profiles to the retrieved values for T and R_eq from our slab retrievals in Fig. 3. As stated in Sect. 2.2, the slab models retrieve an emitting area A that can be converted into an equivalent emitting radius assuming $A=\pi R_{\mathrm{e}\mathrm{q}}^2$ though the emission can also originate from an annulus further out in the disk instead. This can have a significant effect on where the data point will end up on the x-axis of Fig. 3. At a minimum, the inner starting radius of the DALI models, which was set to the dust sublimation radius at 0.08 au, could set the inner radius of such an annulus. However, as can be seen from the CO₂ abundance map presented in Fig. A.2, this may not be sufficient in all cases, as the CO₂ emission mainly originates from further out in the disk where it becomes more abundant. Hence, we parameterized the retrieved emitting areas from our slab retrievals as an annulus starting at the inner edge of the corresponding 70% CBF. For the single-temperature H₂O fits, we matched the fit done at the shortest wavelengths to the CBF of the highest E_up line, and the subsequent wavelength ranges to the lower E_up lines. For the three-temperature-component MCMC, the three components were parameterized as three annuli starting at the inner edge of the CBF of the highest E_up line.

As can be seen from Fig. 3, the agreement between the retrieved temperature from the slab fits and the temperature in the emitting layer of the DALI model is generally quite good, for both the full chemistry and parameterized model. It is clear that the three single-temperature fits (indicated with crosses in Fig. 3) in the 10–14, 13.5–17.5, and 21–24 μm regions mainly trace the “warm,” ~500 K component of the emission (though the 21–24 μm fit in the parameterized model better reflects the cold temperature of this emission, likely due to it being quite strong in this model). While the different wavelength regions do show a slight difference in temperature as expected (colored crosses), they do not capture the full range of temperatures contributing to the emission. Instead, the three-temperature-component MCMC fit does a better job at capturing the gradient in temperature, especially toward the lower temperatures (white crosses). This indicates that a single temperature is not enough to capture the full complexity of the H₂O spectrum even in narrower wavelength ranges, due to the wide range of E_up and A_ij that are present across the rotational spectrum (see, e.g., Meijerink et al. 2009; Banzatti et al. 2025). Thus, H₂O retrievals should be done using at least a multi-temperature fit. Closer inspection of the retrieved slab temperatures shows that, while they generally agree well with the temperature in the emitting layer, they are slightly colder. This could indicate that the lines are slightly sub-thermally excited (see, e.g., Meijerink et al. 2009), as the excitation in the DALI model is calculated in non-LTE whereas the retrieval assumes LTE. We explore the implications of LTE versus non-LTE excitation further in Appendix A.2, and find that assuming LTE excitation in our DALI model indeed leads to a slightly higher retrieved temperature, though the difference is relatively small (~100 K).

Fig. 3 Temperature as a function of radius within the 70% emitting region of the H2O 177,10−164,13 (6371 K), 113,9−100,10 (2438 K), and 83,6−70,7 (1447 K) lines (solid red, pink, and blue lines) for the model using the chemical network (left panel) and the model with parameterized abundances (right panel). The shaded regions represent the minimum and maximum temperature within the emitting region. The red, pink, and blue crosses represent the retrieved T and Req from the single-temperature slab fits in the 10–14, 13.5–17.5, and 21–24 μm regions. The white plus symbols represent the retrieved T and Req values from the three-temperature-component MCMC routine for H2O.

3.2.2 Column density

Fig. 4 compares the retrieved column densities from the LTE slab models to the vertically integrated column density. The retrieved values indicate that the emission traces close to the dust τ = 1 surface for our fiducial models. In Appendix A.3, we discuss how these results differ for a set of models with a higher gas-to-dust ratio. In this case, the dust τ = 1 surface is located deeper in the disk, whereas the emission does not seem to trace much deeper, given the column densities that are retrieved (see Fig. A.7). This indicates that, in cases of a high gas-to-dust ratio, the lines can already become optically thick before the dust does, and thus may not trace the dust τ = 1 surface in that case.

3.3 Investigating the cold H₂O component

Fig. 2 demonstrates that the cold H₂O emission in the model with parameterized abundances is very strong compared to the warmer emission, whereas that is not the case for the model with the full chemistry. We explore here why this could be the case.

3.3.1 Influence of the abundance structure

The only difference between the full chemistry and parameterized abundance models is the abundance structure. Comparing the left panels of Fig. 1, it is clear that the parameterized model has a much higher H₂O abundance above the snow surface with our choice of n(H₂O)/n_H = 10⁻⁴. The lower H₂O abundance in the full chemistry model is due to the efficient photodissociation of H₂O in this layer, as the temperature is only high enough for efficient gas-phase formation (>300 K) in a small sliver at the surface. Thus, the H₂O column this far out in the disk cannot efficiently self-shield against destruction. Both models adopt the same temperature structure, so this cannot be the cause of the observed differences. In fact, a parameterized model with the assumption of T_gas = T_dust (as typically done in previous studies; e.g., Bruderer et al. 2015; Bosman et al. 2017) reveals only a negligible change to the strength of the cold emission (an increase of ~10%), whereas the warm and hot emission is affected much more strongly.

We therefore further investigated the impact of the abundance structure on the H₂O emission. We ran two additional parameterized models where the same abundance structure was adopted as before, but now the H₂O abundance beyond 1 au was lowered by 1 or 2 orders of magnitude to 10⁻⁵ and 10⁻⁶, which allowed the models to more closely match the abundance structure seen in the full chemistry model. The H₂O abundance structure of these new models, along with the fiducial full chemistry and parameterized models, are presented in Fig. A.1. The contours represent the 70% CBFs of the three aforementioned lines. The emitting regions of the two higher-E_up lines remain relatively unchanged between the three parameterized models, but the emitting region of the coldest line is clearly affected by the lower abundance in the outer surface layers. The emitting region of the warm H₂O (pink contours in Figs. 1 and A.1) also increases between the parameterized and full chemistry models, as the parameterized structure now allows it to be abundant at slightly lower temperatures than those at which it would otherwise form efficiently.

Fig. 5 clearly demonstrates the effect of lowering the H₂O abundance above the snow surface by showing four lines of interest between 23.7 and 24 μm, which are a key diagnostic for cold H₂O. Two of these lines (the 9_8,1−8_7,2 and 11_5,6−10_4,7 lines) have an E_up of ~3000 K, whereas the other two lines (the 8_3,6−7_0,7 and 8_4,5−7_1,6 lines) have a lower E_up of ~1500 K. We note that the emitting region of the 8_3,6−7_0,7 line is shown in Fig. 1 in dark blue, and the other 1500 K line traces a nearly identical region in the disk. The two 3000 K lines trace a very similar region as the 2438 K line shown in pink in Fig. 1, though their emission is concentrated slightly higher up (similar to the vertical extent of the 1500 K lines). The flux ratio of these four lines is very sensitive to the temperature of the gas (see, e.g., Banzatti et al. 2023, 2025; Temmink et al. 2024), with the 1500 K lines increasing rapidly in strength with respect to the 3000 K lines when the temperature of the gas is decreased. Additionally, Banzatti et al. (2025) point out that an asymmetry between the 8_3,6−7_0,7 and 8_4,5−7_1,6 lines can constrain the temperature of the coldest emission even further, with a line ratio of ≳1.2 requiring gas temperatures ≲220 K.

We normalized the synthetic spectra of all four models to the 9_8,1−8_7,2 line between the two 1500 K lines. In the full chemistry model (shown in black), all four lines have a similar flux, whereas in the fiducial parameterized model (shown in blue), the colder ~1500 K lines are stronger and show a prominent asymmetry. The strength of these cold lines with respect to the warmer 3000 K lines is reduced as the H₂O abundance above the snow surface is lowered (pink and red lines). We note that these findings do not change when LTE instead of non-LTE excitation is considered. As such, it seems that a relatively high H₂O abundance above the snow surface (≳10⁻⁵) beyond ~1 au is required for the spectrum to show prominent cold H₂O emission and be consistent with several observations (see Banzatti et al. 2023, 2025; Temmink et al. 2024; Romero-Mirza et al. 2024).

In Appendix A.4, we demonstrate that such a high H₂O abundance above the snow surface at large radii (>1 au) is very hard to attain through gas-phase formation alone. The temperature in this layer is simply too low, as is also reflected by the fact that such a reservoir is not present in our fiducial chemistry model (see Fig. 1). Fig. A.10 demonstrates that a strong increase in gas temperature in this layer can produce abundant H₂O above the snow surface (mimicking the abundance structure in our parameterized model); however, this does not come paired with a strong enhancement of the cold H₂O with respect to the hot and warm emission. Hence, this implies that the prominent cold H₂O seen in observations requires a high H₂O abundance above the snow surface combined with a low temperature ≲300 K.

Fig. 4 Vertically integrated H2O column density as a function of radius. The solid lines show the total model column density, and the dashed lines show the model column density integrated up to the dust τ = 1 surface at 15 μm. The red, pink, and blue crosses represent the retrieved N and Req from the single-temperature slab fits in the 10–14, 13.5–17.5, and 21–24 μm regions. The plus symbols in the left column represent the retrieved N and Req values from the three-temperature-component MCMC routine. The left panel depicts the model using the chemical network and the right panel depicts the model with parameterized abundances.

Fig. 5 Synthetic H2O spectra of the fiducial model using the chemical network (black) and three models with parameterized abundances (blue, pink, and red) between 23.7 and 24 μm. The flux is normalized to the 98,1−87,2 line at 23.87 μm.

3.3.2 Temperature diagnostics

To gain further, more comprehensive insight into the temperature distribution probed by the H₂O emission from our models, we created a diagnostic diagram following the method described in Banzatti et al. (2025). In this work, the authors define several clean, non-blended diagnostic H₂O lines with E_up of 6000, 3600, and 1500 K. Their line ratios can then provide useful insight into the relative strengths of different temperature components in the H₂O spectrum. The 3600/6000 K line ratio captures the relative strength of the warm (~400 K) component with respect to the hot (~850 K), and the 1500/3600 K provides the relative strength of the colder component (~200 K) with respect to the warm component. Plotting these two ratios against each other can then, at a glance, demonstrate what temperature component a spectrum (be it a model or an observation) is dominated by.

The temperature diagnostic diagram is presented in Fig. 6, where the samples presented in Banzatti et al. (2025) and Temmink et al. (2025) are shown in gray plus symbols and blue crosses, respectively. EX Lup, as presented in Smith et al. (2025), is indicated using pink triangles. The parameterized abundance models are shown in blue diamonds and the chemistry models are shown in red stars. For both sets of models, the fiducial model is indicated as the largest marker with a thicker outline. Models with a higher fraction of large grains (f_ℓ) and thus a higher gas-to-dust ratio in the IR emitting region (see also Appendix A.3) are also shown. These models lie to the lower-left compared to the fiducial model, indicating that they have a stronger warm component and a weaker cold component compared to the fiducial models. This can be explained by the increased temperature in the surface layers due to the increased dust settling in these models.

The two additional parameterized models presented in Fig. A.1 lie directly below the fiducial parameterized model in Fig. 6, indicating that they have a weaker cold component, but their warm and hot components are unaffected. This demonstrates that the strength of the cold H₂O component in a spectrum is directly linked to the H₂O abundance above the snow surface. The model with increased gas temperature (see Fig. A.10) shows only slightly colder H₂O compared to the fiducial chemistry model, and does not come close to matching the parameterized models, even though the abundance structure is similar. We note that all models lie further to the right compared to most of the data. This could be interpreted as evidence for a high gas-to-dust ratio in the upper layers of these disks. However, it is also possible that our models lack sufficient hot (≳800 K) H₂O. Banzatti et al. (2025) demonstrate the presence of high E_up lines (>7200 K) throughout the rotational H₂O spectrum, which are not present in our synthetic spectra due to the E_up cutoff at 7200 K in the molecular data used due to the lack of collisional rate coefficients for these levels (Faure & Josselin 2008). The inclusion of these lines may enhance the hot component and allow our models to better match the observations.

In Fig. 6, several disks that are known to exhibit a strong cold H₂O component in their spectra are labeled: IRAS-04385 (Banzatti et al. 2025), FT Tau, and XX Cha (Temmink et al. 2025). The 1500/3600 K line ratio of FT Tau and XX Cha – two disks, which are noted to have strongly enhanced cold H₂O – generally agrees with our parameterized models, thus potentially indicating that these sources indeed have a relatively high H₂O abundance above the snow surface. The same holds true for EX Lup, which was also demonstrated to have strong cold H₂O by Smith et al. (2025), who speculate that this may be the result of an outburst (see also Sect. 4.3).

Comparing our models to IRAS-04385 brings up another interesting finding: its 1500/3600 and 3600/6000 K line ratios match our fiducial full chemistry model very closely. However, their spectra are visually quite different, especially at 23.8 μm. In IRAS-04385, the two 1500 K lines are much stronger than the 3000 K line (resembling our fiducial parameterized model, the blue line in Fig. 5; see Fig. 11 in Banzatti et al. 2025), whereas this line ratio is relatively flat in our fiducial full chemistry model (black line in Fig. 5). This indicates that the 1500/3600 K line ratio is not always directly indicative of what the lines at 23.8 μm look like, and some caution should be taken before it is definitively concluded whether a source has strong cold H₂O emission or not.

Fig. 6 H2O temperature diagnostic diagram, as defined in Banzatti et al. (2025). The observations presented in Banzatti et al. (2025), Temmink et al. (2025), and Smith et al. (2025) are shown as gray plus symbols, blue crosses, and pink triangles, respectively. The parameterized models (denoted as PA) presented in this work are shown as blue diamonds, and the chemistry models are shown as red stars. In both cases, the fiducial model is indicated as the largest symbol with a thicker black outline.

Fig. 7 H2O abundance maps of four models using the chemical network with varying accretion luminosity. The τ = 1 surface of the dust continuum at 15 μm is indicated with a dashed blue line. Red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371K), 113,9−100,10 (Eup = 2438K), and 83,6−70,7 (Eup = 1447K) lines, respectively.

Fig. 8 Left: Input blackbody spectra for an L* = 0.8L⊙ source with an accretion luminosity of Lacc = 0, 0.002, 0.02 and 0.2L⊙ (black, yellow, pink, and blue). Right: Synthetic H2O spectra corresponding to the four models in the left panel.

4 Discussion

The analysis presented in Figs. 5 and 6 indicates that a high H₂O abundance (≳10⁻⁵) above the snow surface at radii ≳1 au is required for prominent cold H₂O emission to be present in the spectrum. However, such a high abundance is not naturally present in our full chemistry model as the temperature is too low and H₂O is efficiently photodissociated in this region. We explore several possibilities how one can obtain a higher H₂O abundance above the snow surface, to understand how this prominent cold H₂O emission could arise in observations.

4.1 Lower FUV luminosity

First, we investigated whether a lower photodissociation rate could enhance the H₂O abundance above the snow surface, leading to stronger cold H₂O emission. This could occur, for example, due to a lower FUV flux coming from the central source. To this end, a small grid of four models with full chemistry and varying accretion (FUV) luminosity was run. Instead of the AS 209 input stellar spectrum, we used a 4000 K blackbody with an additional 10 000 K component scaled by the mass accretion rate (see Kama et al. 2016; Visser et al. 2018). All other parameters related to the disk structure and input elemental abundances were kept constant. This is summarized in Table A.2. The resulting H₂O abundance structure is presented in Fig. 7. Clearly, the H₂O abundance above the snow surface increases, as the input FUV luminosity of the model decreases due to the decreased photodissociation rate.

The resulting synthetic spectra from the four models are presented in Fig. 8, once again showing the 23.7–24 μm region. As expected, due to the increase in H₂O abundance with decreasing FUV luminosity, the cold 1500 K lines (8_3,6−7_0,7 and 8_4,5−7_1,6) indeed increase in strength with respect to the 3000 K lines (9_8,1−8_7,2 and 11_5,6−10_4,7) with decreasing FUV luminosity. However, the difference in line ratio is not nearly as apparent as in Fig. 5; the H₂O abundance above the snow surface, while it does increase when the input FUV luminosity is lowered, does not exceed n(H₂O)/n_H ~ 10⁻⁶. Thus, the variety in this line ratio as seen in MIRI observations (see, e.g., Banzatti et al. 2023, 2025; Romero-Mirza et al. 2024) cannot be explained solely by a difference in FUV luminosity.

4.2 Vertical mixing and enhanced H₂O self-shielding

Second, we should consider that our thermochemical models are static and thus do not include dynamical processes such as dust transport. Recent work has provided growing evidence that dust dynamics, specifically radial drift, plays an important role in disks that show strong cold H₂O emission. These disks, with strong cold H₂O, are found to be generally more compact in millimeter dust emission than the disks in which the cold H₂O emission is less prominent (Banzatti et al. 2020, 2023, 2025; Romero-Mirza et al. 2024), although whether this generally holds is challenged by new data. Temmink et al. (2025) find that not all of the compact disks in their sample show an enhancement in the cold emission. It is postulated that the cold H₂O emission is a signature of strong radial drift in these disks that brings a large volume of ices across the H₂O snow line. The subsequent sublimation increases the amount of cold vapor mass at the snow line, which translates into an enhanced flux of the coldest H₂O lines. However, our work demonstrates that it is not the reservoir within the midplane snow line that the coldest H₂O emission is most sensitive to. Instead, it is the reservoir above the snow surface at larger radii that impacts these lines the most. Hence, there needs to be a contribution from vertical mixing that carries the icy grains upward across the 2D snow surface, in conjunction with the radial drift.

For vertical mixing to enhance the H₂O abundance to the degree that seems to be required from our modeling (reaching a fractional abundance of ≳10⁻⁵), it needs to overcome the fast photodissociation rates in the layer above the snow surface. In the left panel of Fig. 9, we present the H₂O photodissociation timescale for the fiducial chemistry model. In the right panel, we present a simple calculation of the turbulent vertical mixing timescale, following Xie et al. (1995) and Semenov et al. (2006). These works assume the vertical mixing to be set by turbulent diffusion, thus allowing the mixing timescale to be approximated as t_mix ≈ H²/D_turb. Here, H = c_s/Ω_K is the scale height of the disk and $D_{\text{turb\hspace{0.17em}}}=\alpha c_{\mathrm{s}}^2/{\mathrm{\Omega }}_{\mathrm{K}}$ is the turbulent diffusion coefficient. c_s denotes the sound speed (assumed to be isothermal), Ω_K denotes the Keplerian angular frequency and α denotes the viscosity strength. We assumed a strong turbulence of α = 10⁻² for our calculation. We stress that this is a best-case scenario, as many observations have estimated the turbulence in disks to be much lower (≲10⁻³; see, e.g., Flaherty et al. 2020; Villenave et al. 2025), with only a few exceptions (e.g., Paneque-Carreño et al. 2024; Flaherty et al. 2024).

Fig. 9 shows that the photodissociation timescale is short (≲1 yr; seen also in, e.g., Xu et al. 2019; Kanwar et al. 2024) in the surface layer above the snow surface where the H₂O abundance is correlated with the cold H₂O emission. The mixing timescale, on the other hand, is much longer (10–100 years), even in the best-case scenario with very strong turbulence. Thus, this simple timescale comparison does not make it seem likely that vertical mixing by turbulent diffusion will be able to enhance the cold H₂O reservoir above the snow surface by a significant amount.

Still, a comparison of timescales alone may not provide the full picture. The effects of vertical mixing by turbulent diffusion have been implemented in the thermochemical code ProDiMo by Woitke et al. (2022) by iteratively solving the reaction-diffusion equations until convergence is reached. They show that this enhances the H₂O abundance in the uppermost surface layers (T ≳ 1000 K) of the disk by up to 1–2 orders of magnitude (though their models do not reach fractional abundances higher than ~10⁻⁶ at 1 au or beyond). These models do not consider H₂O UV-shielding (Bethell & Bergin 2009), radial mixing, or the effects of dust transport, so the continuous replenishment of ices from the outer disk could perhaps enhance these effects even further.

Moreover, the simple timescale comparison does not take into account the fact that dynamical processes such as radial drift or vertical mixing may build up a sufficient column of H₂O for UV-shielding to become important. This process was accounted for during the calculation of the chemistry and thermal balance in DALI (Bosman et al. 2022a,b); however, DALI cannot account for the extra H₂O vapor brought to the inner disk by dynamical processes. Therefore, in reality, the effects of H₂O UV-shielding may be more efficient than what is accounted for in our models, further reducing the photodissociation rate and thus allowing for a more significant buildup of cold H₂O vapor above the snow surface.

A simple estimate reveals that reducing the photodissociation timescale by a factor ~100 by self-shielding (which would bring it a lot closer to the mixing timescale) would require a H₂O column density of ~10¹⁸ cm⁻² (taking an average photodissociation cross section for H₂O of σ_av ~ 5 × 10⁻¹⁸ cm²; Fillion et al. 2004). This value is not unreasonable given the constraints obtained from recent observations. Temmink et al. (2024) fit the rotational H₂O spectrum of DR Tau (where the 1500 K lines at 23.8 μm are quite strong), retrieving a column density of the coldest component between log N = 17.6–17.9, depending on their exact method. The column density profiles derived by Romero-Mirza et al. (2024) also generally come close to this value for the regions in the disk where T ≲ 300 K.

Thus, it seems that a clear consensus on whether vertical mixing can produce the H₂O abundance needed to produce strong cold H₂O emission would require more sophisticated dynamical modeling. The timescales by themselves suggest that the photodissociation in the disk’s surface layers is too fast to allow for a significant buildup of H₂O vapor by vertical mixing; however, a more thorough exploration of disk dynamics and its effects on the IR H₂O spectrum is warranted for future work. Greenwood et al. (2019) has shown that dust evolution by itself will enhance IR emission over time due to a decreased opacity from small micron-sized grains in the surface layers, though they do not include the delivery of volatiles to the inner disk. Alternatively, one could also consider a mechanism, such as surface accretion flows, that can bring small icy grains inward across the snow surface via the surface layers, rather than being mixed upward from the midplane (e.g., Bitsch et al. 2014).

Fig. 9 H2O photodissociation timescale and turbulent diffusion mixing timescale of the fiducial model using the chemical network. The 1,10 and 100 year timescale contours are indicated in white in both panels. The dark blue contours in both panels indicate the 10−5, 10−6, and 10−8 H2O abundance contours of the fiducial chemistry model (see also top left panel of Fig. 1). The light blue line represents the H2O snow surface.

4.3 Accretion outbursts

Finally, Smith et al. (2025) provide a promising alternative explanation that links the strength of cold H₂O and the H₂O abundance above the snow surface: accretion outbursts. The authors demonstrate that EX Lup, an M0 star with known accretion variability, experienced a cold H₂O vapor burst during its accretion outburst in 2008. This is observed by an increase in the strength of the cold, low-energy H₂O lines between two Spitzer epochs, which the authors attribute to the recession of the snow line to larger radii, allowing a large volume of ice to sublimate. A comparison with more recent JWST/MIRI observations shows that this strong cold emission is still present, suggesting long (>10 yr) refreeze-out timescales.

Appendix D of Smith et al. (2025) demonstrates using ProDiMo thermochemical models how the H₂O abundance structure can be altered by an accretion outburst. The recession of the snow line creates a layer of highly abundant H₂O (n(H₂O)/n_H ~ 10⁻⁴) right at the snow surface. This coincides very well with our predictions in Sect. 3.3, where we demonstrate that an increase in H₂O abundance in this layer is directly responsible for an increase in the strength of the cold H₂O lines. Moreover, their Figure 10 demonstrates that the outburst increases the H₂O column density at large radii by up to 4 orders of magnitude, reaching values of ~10²¹ cm⁻², meaning that this reservoir would be able to self-shield against photodissociation quite effectively.

We demonstrate a similar finding in Fig. A.10, where a strong increase in gas temperature (perhaps due to an outburst) indeed leads to an enhancement in the H₂O abundance above the snow surface, resembling our fiducial parameterized model and the model presented in Smith et al. (2025). In our case, this is purely caused by the enhancement in gas-phase H₂O formation due to the increased temperature. In our own model grid with varying accretion luminosity (Sect. 4.1), we assumed the accretion state to be static and thus evolved the chemistry for 1 Myr (during which the surface layers already reach a steady-state chemistry). Hence, we did not find the accretion luminosity to have such a drastic effect on the H₂O abundance on longer timescales. In contrast, Smith et al. (2025) evolved their outburst model, starting from a quiescent state, for only a very short amount of time (half a year), causing the strong change in the H₂O abundance structure. This demonstrates that processes occurring on short timescales may also leave an imprint in the IR H₂O spectrum.

Appendix A.4 further demonstrates that the low temperature in the layer above the snow surface in the fiducial chemistry model does not allow H₂O to form abundantly (i.e., reaching n(H₂O)/n_H ~ 10⁻⁴; see Fig. A.9). This only becomes more feasible if the temperature is increased (e.g., due to an accretion outburst). However, since the latter also enhances the emission from the hot and warm H₂O, this does not lead to a spectrum dominated by cold H₂O. Instead, this implies one of two conclusions: 1) the H₂O needs to remain abundant in this layer while the disk cools down again (potentially aided by UV shielding), or 2) the H₂O needs to become abundant through a mechanism other than gas-phase formation, such as H₂O delivery from radial drift followed by vertical mixing. Given the rapid photodissociation timescales in the surface, the latter case possibly implies that some sort of continuous replenishment is required.

Thus, accretion outbursts seem to produce a temporary burst of cold H₂O vapor that can be traced in the IR with, for example, JWST/MIRI. However, it remains debated whether this explanation is generally applicable to all disks that exhibit cold H₂O, since not all disks show strong evidence of accretion variability. Moreover, it seems that a rather large outburst is required, as Smith et al. (2025) did not find the moderate outburst of EX Lup in 2023 to cause a significant change in the spectrum.

Fig. 10 Synthetic H2O spectra of the fiducial model using the chemical network (black) and the fiducial parameterized model (abundance above the snow surface of 10−4; blue) between 23 and 33, 40 and 55, 60 and 80, and 120 and 185 μm. The flux is normalized to the 98,1−87,2 line (2892 K) at 23.87 μm. Blue triangles indicate H2O lines with Eup < 2500 K.

4.4 Future outlook in the far-IR

Current observations and theoretical efforts focus on the IR H₂O spectrum of disks as seen with JWST/MIRI and how dynamical processes may affect this (e.g., Banzatti et al. 2023, 2025; Romero-Mirza et al. 2024; Kalyaan et al. 2021, 2023; Mah et al. 2023). While JWST has been able to provide new insight into the warm (~500 K) H₂O content of disks, lower E_up lines tracing the colder content are only scarcely available in the spectrum, with no lines with E_up < 800 K present in the MIRI wavelength range. Additionally, these lines are located mostly in the longest wavelength channel, which has the lowest sensitivity. To fully understand the cold H₂O content in disks, the material that is transported across the snow line by radial drift, observations at longer wavelengths, in the far-IR, are crucial (see also, e.g., Zhang et al. 2013; Pontoppidan et al. 2018; Kamp et al. 2021).

While the far-IR window is currently unavailable for disk spectroscopy, the future looks promising, with both a confirmed balloon mission and a space mission under development that will be able to survey the cold H₂O content in planet-forming disks. The Planetary Origins and Evolution Multispectral Monochromator (POEMM) stratospheric balloon mission is a 2 m telescope with a high-resolution (R = 10⁵) spectrograph covering wavelengths from 35–112 μm (Stacey et al. 2024). One of its main science goals is to characterize the H₂O emission from Herbig and T Tauri disks, to further our understanding of planet formation. The PRobe far-Infrared Mission for Astrophysics (PRIMA) mission is a 1.8 m space telescope, with the FIRESS instrument providing spectroscopy of the 24–235 μm range at a superior spectral resolution to JWST/MIRI (R = 4400 × 112 μm/λ), and with a better sensitivity than POEMM. It builds technically and scientifically on earlier design concepts such as the SPace Infrared telescope for Cosmology and Astrophysics (SPICA) mission (Roelfsema et al. 2018; Kamp et al. 2021).

Fig. 10 presents the synthetic spectra for the fiducial chemistry and parameterized models (black and blue lines). The fluxes are normalized to the flux of the 9_8,1−8_7,2 line at 23.87 μm as also done in Fig. 5 but are now shown from 23 to 185 μm. This illustrates that there are many more lines in the far-IR that are very sensitive to the cold H₂O abundance above the snow surface. Hence, these lines may also make for good tracers of the amount of H₂O ice crossing the snow line. Additionally, the two missions highlighted above will also provide much higher spectral resolution (PRIMA reaches R ~ 20 000 at λ = 25 μm), and will thus be able to provide kinematical information as well for the brightest sources, allowing for the spatial distribution of the emission, and thus the location of the snow line, to be constrained.

5 Conclusions

In this study, we tested several retrieval techniques that are commonly used to interpret JWST observations of H₂O on synthetic DALI spectra. We summarize our findings from these tests as follows:

• Single-temperature slab fits generally trace the warm (~500 K) component of the H₂O emission, even when the fits are performed in different narrow windows within the MIRI wavelength range. This indicates that, even within a relatively narrow window, a single temperature is not enough to capture the full complexity of the H₂O emission;

• A three-component MCMC fit to the entire MIRI wavelength range (following, e.g., Temmink et al. 2024) is better suited to capture the full gradient in temperature present in the H₂O spectrum, as it also captures the contributions from hotter (>600 K) and colder (<300 K) emission;

• Retrieved temperatures generally agree well with the model, though they may slightly underestimate the true temperature in the emitting layers due to non-LTE effects such as sub-thermal excitation;

• The retrieved column density, both from single- and three-component fits, generally traces close to the dust τ = 1 surface for our fiducial models. However, in models with more settled dust where the τ = 1 surface lies deeper into the disk. The lines become optically thick before the dust does, and the emission therefore does not trace all the way down to the τ = 1 surface;

• The retrievals on CO₂ emission yield the same conclusions as above, and we find that the ¹³CO₂ emission retrieves a lower temperature than ¹²CO₂ due to it tracing deeper into the disk.

Additionally, we further investigated the cold (~200 K) component of the H₂O emission to understand the origin of this emission. We summarize our findings as follows:

• We demonstrate that the strength of the cold H₂O emission is directly linked to the H₂O abundance above the snow surface at larger radii (≳1 au). This implies that the disks in which strong cold H₂O has been observed (e.g., Banzatti et al. 2023, 2025) likely have a high abundance (≳10⁻⁵) in this reservoir – higher than that predicted by our fiducial full chemistry model;

• We demonstrate that the temperature needs to be low (<300 K), meaning this high abundance must be achieved through a mechanism other than gas-phase formation, which is not efficient at these temperatures in gas exposed to intense UV radiation;

• A lower FUV luminosity from the central star, for example, due to a low mass accretion rate, can increase the H₂O abundance above the snow surface. However, this effect is too small to explain the diversity in cold H₂O strength seen in observations;

• Instead, vertical mixing of ices across the snow surface, in conjunction with their inward transport by radial drift, could possibly enhance the H₂O abundance beyond what is considered in our static thermochemical models. This would fall in line with recent work tentatively identifying cold H₂O emission as a signature of strong radial drift (Banzatti et al. 2023, 2025; Romero-Mirza et al. 2024), though further investigation is needed to confirm this;

• Accretion outbursts provide an interesting alternative explanation for enhanced cold H₂O, as Smith et al. (2025) demonstrate. These outbursts may enhance the H₂O abundance above the snow surface, which, as our study has demonstrated, directly correlates with the strength of cold emission.

With this work, we provide further indications from a modeling perspective that the cold H₂O emission seen with JWST may be a signature of a dynamical process, as argued previously from an observational perspective. To understand the origins of the cold H₂O emission, and to give us further insight into how JWST observations could possibly probe the H₂O at its snow line, the effects of dust transport and mixing on the observed IR emission should be investigated in future work.

Acknowledgements

We thank the referee for their constructive feedback that helped improve this manuscript. Astrochemistry in Leiden is supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101019751 MOLDISK). E.v.D. also acknowledges support from the Danish National Research Foundation through the Center of Excellence “InterCat” (DNRF150).

Appendix A Additional model explorations

A.1 Retrievals on CO₂ emission

Since CO₂ emission is seen in many disks observed with JWST/MIRI (e.g., Grant et al. 2023; Vlasblom et al. 2025; Salyk et al. 2025), we also investigate how the temperature and column density retrieved from LTE slab fits on synthetic CO₂ spectra compare to our models. Just as for H₂O, we perform our retrievals both on the model using the full chemistry, and on a model with a parameterized abundance. The abundance maps for these models are shown in Fig. A.2. The 70% emitting regions of the ¹²CO₂ and ¹³ CO₂ ν = 1–0 Q(20) lines (green and purple, respectively) are also indicated. These species trace further out into the disk than the H₂O lines (11_3,9−10_0,10 line with E_up = 2438 K is indicated in pink), and the ¹³CO₂ traces slightly deeper into the disk than the ¹²CO₂, since it is optically thinner. The reason that CO₂ traces further out than H₂O is most likely due to CO₂ being more efficiently formed at lower temperatures, where the OH + CO ⟶ CO₂ + H reaction is more efficient than the OH + H₂ ⟶ H₂O + H reaction (Charnley 1997; van Dishoeck et al. 2013; Walsh et al. 2015).

Additionally, the version of DALI used in this work includes the effects of UV shielding by H₂O (Bethell & Bergin 2009) in its calculation of the chemistry and thermal balance, as implemented in Bosman et al. (2022a). This process is most efficient in the innermost regions of the disk, where the H₂O column density can build up substantially due to efficient gas-phase formation. As a result, the formation of CO₂ is strongly quenched in these regions, as the self-shielding of H₂O prevents the production of OH via H₂O photodissociation, thus preventing CO₂ production. Additionally, the photodissociation rate of CO₂ is higher in these innermost regions as it cannot benefit from self-shielding at these column densities. This creates a gap in the CO₂ abundance within 0.4 au, as first demonstrated in Bosman et al. (2022b), and also seen in our model (left panel of Fig. 1).

For consistency, this effect is recreated in the parameterized model; thereby, we use a more complex abundance structure to better match the model with the full chemistry. We define the high-abundance reservoir as the region where 200 K > T_dust > 70 K and A_V > 5 mag, where we set the fractional abundance to n(CO₂)/n_H = 5 × 10⁻⁵. We also define a lower-abundance surface layer as the region where 1 < A_V < 5 and T_dust > 70 K, where we set the fractional abundance to n(CO₂)/n_H = 10⁻⁷. The fractional CO₂ abundance in the rest of the disk is set to 10⁻¹⁰.

Since ¹²CO₂ and ¹³CO₂ each only have one main emission feature in the MIRI wavelength range, we fit these features with a single LTE slab model between 13.5 and 17.5 μm, retrieving a single best-fit value for T, N and R_eq. The best-fit model was obtained from a grid of models using the χ² method described in Sect. 2. In this grid, T was varied linearly between 100 and 1000 K in steps of ΔT = 22.5 K and N was varied in log-space between 10¹⁴ and 10²³ cm⁻² with steps of Δ log N = 0.225. For ¹³CO₂, T was varied between 50 and 500 K in steps of ΔT = 11.25 K.

The results of these fits are shown in Fig. A.3, and just as for H₂O, the agreement between the retrieved temperature and the temperature within the emitting layer of the DALI model is quite good. We find that the ¹³CO₂ emission traces a slightly lower temperature than the ¹²CO₂, likely due to it tracing deeper down into the disk. The retrieved column densities are presented in Fig. A.4, which trace close to the τ = 1 surface, just as was found for H₂O.

A.2 LTE versus non-LTE excitation

As mentioned in Sect. 2, the excitation in DALI can be calculated either in LTE or non-LTE, if collisional data are available for the species in question. In the main body of this work, the excitation for H₂O is calculated in non-LTE, from which the synthetic spectra are then generated. However, all subsequent retrievals on these synthetic spectra assume LTE excitation. To understand the implications of this, we calculate the excitation for our two fiducial models (full chemistry and parameterized abundances) in LTE as well, generate new synthetic spectra and run our LTE retrievals again. The results are presented in Fig. A.5.

Fig. A.5 shows the retrieved temperatures of our slab fits in three wavelength ranges (10-14, 13.5-17.5, and 21-24 μm) for both our LTE (opaque crosses) and non-LTE (faded crosses) excitation models. The temperature retrieved from our LTE model is slightly higher than the temperature retrieved from our non-LTE (fiducial) model, both in the full chemistry and parameterized case. We note that this effect is also seen for the three-component MCMC fit, though this is not shown on Fig. A.5. This clearly indicates that, if non-LTE effects are accounted for, the H₂O emission seen with JWST is slightly sub-thermally excited, as confirmed in previous work (Meijerink et al. 2009). Hence, the retrieved temperature from an observation can in reality be slightly too low (or the emitting radius should be slightly larger), though the discrepancy is relatively small at only ~100 K.

A.3 Increased gas-to-dust ratio

It is predicted that the gas-to-dust ratio in the IR emitting layer is elevated above the canonical ISM value of 100 (Meijerink et al. 2009). In DALI, this can be parameterized either by the gas-to-dust mass ratio (which is kept at 100 for all models in this work), or by the parameter f_ℓ, which sets the mass fraction of the large grain population. Since the large grains have a reduced scale height with respect to the small grains (which follow the vertical extent of the gas), a larger value for f_ℓ increases the gas-to-dust ratio in the upper layers of the disk. In the models presented in the main body of this work, f_ℓ is set to 0.9, which yields a gas-to-dust ratio of ~10³ in the IR emitting layers. We run both fiducial models (full chemistry and parameterized abundances) with f_ℓ = 0.99 and 0.999, yielding even higher gas-to-dust ratios of 10⁴ and 10⁵ in the upper layers of the disk.

Fig. A.6 presents the retrieved temperatures for the f_ℓ = 0.999 models. They largely follow the same conclusions as the fiducial models: the single-temperature fits mainly trace the warm ~500 K component, whereas the three-component MCMC allows for the full temperature gradient present in the emitting region to be captured better. We note that the full chemistry model with f_ℓ = 0.999 is quite warm, and thus the contribution from the coldest (<300 K) component is very weak, meaning that this component of the MCMC is poorly constrained.

Fig. A.7 also presents the retrieved column densities. Whereas these mainly traced the τ = 1 surface in the fiducial models, in this model we see that this is not the case. Since the large dust grains are so strongly settled toward the midplane in this model, the τ = 1 surface lays very deep into the disk, and thus the total observable column of gas is large. Consequently, the H₂O gas becomes optically thick before the dust does, and the retrieved column densities no longer trace the dust τ = 1 surface, but rather the gas τ = 1 surface. This is similar to that observed for HCN by Bruderer et al. (2015).

A.4 Influence of the gas temperature

In the main body of this work, we keep the temperature structure of the disk constant between our full chemistry and parameterized models, which allows us to isolate the effects of changes to the abundance structure. Here, we explore how the gas temperature, as well as the strength of the radiation field, influences the H₂O abundance in the IR emitting layer.

For this analysis, we construct a very simple chemical network following Kristensen et al. (2017). We assume that all H₂O is formed through gas-phase reactions with H₂. This is governed by two reactions: O + H₂ ⟶ OH + H and OH + H₂ ⟶ H₂O + H. We also assume that H₂O is only destroyed through photodissociation, where it is rapidly dissociated all the way back to atomic O (this can happen either in a single step or in two steps, where OH dissociates into O + H right after the H₂O photodissociation). The rate of photodissociation is given by R_pd = G₀k_pdn_H2O, where G₀ is the radiation field strength in units of the interstellar radiation field, k_pd = 8.0 × 10⁻¹⁰ s⁻¹ (van Dishoeck et al. 2006) and n_H2O is the H₂O abundance. Thus, assuming equilibrium, we find that $$\frac{\mathrm{d}{n}_{{\mathrm{H}}_2\mathrm{O}}}{\mathrm{d}t}=k_{{\mathrm{H}}_2\mathrm{O}}{n}_{{\mathrm{H}}_2}{n}_{\mathrm{O}\mathrm{H}}-G_0{k}_{\mathrm{p}\mathrm{d}}{n}_{{\mathrm{H}}_2\mathrm{O}}=0.$$(A.1)

Regarding OH, we assume this to be rapidly produced from O + H₂ ⟶ OH + H and destroyed through OH + H₂ ⟶ H₂O + H (we do not consider OH photodissociation), and thus the molecule will not have a substantial abundance with respect to O and H₂O. Thus, in equilibrium, $$\frac{\mathrm{d}{n}_{\mathrm{O}\mathrm{H}}}{\mathrm{d}t}=k_{\mathrm{O}\mathrm{H}}{n}_{{\mathrm{H}}_2}{n}_{\mathrm{O}}-k_{{\mathrm{H}}_2\mathrm{O}}{n}_{{\mathrm{H}}_2}{n}_{\mathrm{O}\mathrm{H}}=0.$$(A.2)

From these two equations, we find that $$n_{{\mathrm{H}}_2\mathrm{O}}=\frac{k_{\mathrm{O}\mathrm{H}}}{k_{\mathrm{p}\mathrm{d}}}\frac{n_{{\mathrm{H}}_2}{n}_{\mathrm{O}}}{G_0}.$$(A.3)

Finally, we assume that n_O = 2.88 × 10⁻⁴ n_H2, following the assumed O/H abundance in our models, and we find that k_OH(T) = 3.13 × 10⁻¹³(T/300 K)^2.7 exp(−3150 K/T) (UMIST22; Millar et al. 2024).

In Fig. A.8, we present the H₂ density, radiation field strength G₀, and gas temperature in our fiducial chemistry model. The solid contour indicates our parameterized H₂O abundance structure, to illustrate what the relevant conditions are in the different parts in the IR emitting layer. For this analysis, we divide the H₂O reservoir into two parts, just as was done in Sect. 3.3: an inner, warm reservoir that is located within 1 au, and an outer, colder reservoir outside 1 au. We only consider the upper, emitting layers above Z/R = 0.15. Fig. A.8 demonstrates that the inner reservoir has a relatively high H₂ density (~10¹⁰−10¹¹ cm⁻³; see van Dishoeck et al. 2021 for a lower density case), a radiation field between 10⁴−10⁶ G₀ and a temperature ≳300 K. The outer reservoir, on the other hand, has a slightly lower H₂ density, a slightly weaker radiation field (between 10³−410⁵ G₀) and a lower temperature ≲400 K.

Fig. A.9 presents the expected H₂O abundance as a function of temperature and G₀ from our calculation, where the red and blue boxes represent the estimated conditions of the inner and outer H₂O reservoir, as specified above. This figure clearly demonstrates that building up abundant H₂O (n(H₂O)/n_H ~ 10⁻⁴) is easily attainable in the inner, warmest regions of the disk (especially since the H₂ density is also higher), but very hard to attain in the outer reservoir, as the temperature is clearly too low.

Thus, Fig. A.10 presents a scenario in which the disk gas temperature has been strongly increased throughout the entire disk, by a factor of 2. We do not increase the dust temperature. Thus, the snow surface remains in the same place. The increase in gas temperature allows H₂O to form abundantly out to much larger radii and the abundance structure closely resembles our parameterized model (see Fig. 1). Still, the spectrum in the third panel of Fig. A.10 demonstrates that the enhancement in the cold H₂O lines is much less than that attained by our parameterized model. After all, while the cold H₂O is slightly enhanced due to an increase in its emitting area, the warm and hot H₂O have undergone a similar increase in their emitting area, and thus the relative strength of the cold lines has not changed much.

Thus, a spectrum dominated by cold H₂O seems to require a high H₂O abundance in the layer above the snow surface, along with a low temperature ≲300 K. Such a reservoir is not present in this model, as the lowest-temperature gas-phase H₂O of the fiducial model at ~150 K is now twice as hot. Since Fig. A.10 has demonstrated that gas-phase formation alone at temperatures below 300 K does not allow abundant H₂O to be built up, this leads to one of two conclusions: 1) the H₂O must remain abundant after the disk cools down from a heating event (such as an accretion outburst; see also Sect. 4.3), which can potentially be aided by UV shielding, or 2) the H₂O must become abundant through some other mechanism than gas-phase formation, such as the sublimation of icy pebbles followed by vertical mixing. The latter may also require some form of continuous replenishment due to the short photodissociation timescales in the upper layers. We note that the effects of additional ice sublimation are not included in this model, as the dust temperature was not enhanced along with the gas temperature. Regardless, the H₂O ice reservoir is located below Z/R = 0.2, below the cold H₂O emitting region and would thus not contribute to its emission unless it is mixed upward.

Fig. A.1 H2O abundance maps of the fiducial model using the chemical network (first panel) and three models with parameterized abundances of 10−4, 10−5, and 10−6 beyond 1 au (second, third, and fourth panels). The τ = 1 surface of the dust continuum at 15 μm is indicated with a dashed blue line. Red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371 K), 113,9−100,10 (Eup = 2438 K), and 83,6−70,7 (Eup = 1447 K) lines, respectively.

Fig. A.2 Abundance maps of CO2 for the fiducial models with fℓ = 0.9. The left panel depicts the model using the chemical network, and the right panel depicts the model with parameterized abundances. The τ = 1 surface of the dust continuum at 15 μm is indicated with a dashed blue line in all panels. The green and purple contours represent the 70% emitting regions of the 12CO2 and 13CO2 ν2 = 1–0 Q(20) lines, respectively, and the pink contours represent the 70% emitting region of the H2O 113,9−100,10 (2438 K) line.

Fig. A.3 Temperature as a function of radius within the 70% emitting region of the 12CO2 and 13CO2 ν2 = 1–0 Q(20) lines (solid green and purple lines). The shaded regions represent the minimum and maximum temperature within the emitting region. The green and purple crosses represent the retrieved T and Req for 12CO2 and 13CO2. The left panel depicts the model using the chemical network, and the right panel depicts the model with parameterized abundances.

Fig. A.4 Vertically integrated CO2 column density as a function of radius. The solid lines show the total model column density, and the dashed lines show the model column density integrated up to the dust τ = 1 surface at 15 μm. The green and purple crosses represent the retrieved N and R for 12CO2 and 13CO2. The left panel depicts the model using the chemical network, and the right panel depicts the model with parameterized abundances.

Fig. A.5 Temperature as a function of radius within the 70% emitting region of the H2O 177,10−164,13 (6371 K), 113,9−100,10 (2438 K), and 83,6−70,7 (1447 K) lines (solid red, pink, and blue lines) for the model using the chemical network (left panel) and the model using parameterized abundances (right panel) with the excitation calculated in LTE. The red, pink, and blue crosses represent the retrieved T and Req from the single-temperature slab fits in the 10-14, 13.5-17.5, and 21-24 μm regions. The opaque crosses represent the fits done on the LTE-excitation models, and the faded crosses represent those done on the non-LTE excitation models (corresponding to the fits shown in Fig. 3).

Fig. A.6 Temperature as a function of radius within the 70% emitting region of the H2O 177,10−164,13 (6371 K), 113,9−100,10 (2438 K), and 83,6−70,7 (1447 K) lines (solid red, pink, and blue lines) for the models with fℓ = 0.999 with the full chemistry (left panel) and with parameterized abundances (right panel). The shaded regions represent the minimum and maximum temperature within the emitting region. The red, pink, and blue crosses represent the retrieved T and Req from the single-temperature slab fits in the 10-14, 13.5-17.5, and 21-24 μm regions. The white plus symbols represent the retrieved T and Req values from the three-temperature-component MCMC routine.

Fig. A.7 Vertically integrated H2O column density as a function of radius for the models with fℓ = 0.999 with the full chemistry (left panel) and with parameterized abundances (right panel). The solid lines show the total model column density, and the dashed lines show the model column density integrated up to the dust τ = 1 surface at 15 μm. The red, pink, and blue crosses represent the retrieved N and Req from the single-temperature slab fits in the 10-14, 13.5-17.5, and 21-24 μm regions. The white plus symbols represent the retrieved N and Req values from the three-temperature-component MCMC routine.

Fig. A.8 Maps of the H2 density, UV radiation field G0, and gas temperature of the fiducial chemistry model. The solid black (and white in the second panel) contours represent the H2O abundance structure of the fiducial parameterized model (see Fig. 1).

Fig. A.9 H2O abundance as a function of gas temperature and UV radiation field, following a simple chemical model (see Sect. A.4). The red and blue boxes represent the conditions within the inner and outer H2O reservoir.

Fig. A.10 H2O abundance, gas temperature, and H2O spectrum between 23.7 and 24 μm for a model with a gas temperature increased by a factor of 2 with respect to the fiducial chemistry model. In the first panel, red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371K), 113,9−100,10 (Eup = 2438K; also shown in black in the second panel), and 83,6−70,7 (Eup = 1447K) lines, respectively. The dust τ = 1 surface at 15 μm is indicated with a dashed blue line in the first and second panel. In the third panel, the fiducial chemistry (chem.) model is shown in black; the model with increased temperature is shown in red; and the fiducial parameterized model (denoted as PA) is shown in blue.

References

All Tables

All Figures

Fig. 1 Abundance maps of H2O (left and middle panels) and temperature map (right panel) for the fiducial models with fℓ = 0.9. The left panel depicts the model using the chemical network, and the middle depicts the model with parameterized abundances. In all panels, the dust τ = 1 surface at 15 μm and the AV = 1.5 surface are indicated with a dashed blue line and a dotted black line, respectively. The red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371 K), 113,9−100,10 (Eup = 2438 K), and 83,6−70,7 (Eup = 1447 K) lines, respectively, representing hot, warm, and cold H2O. The warm H2O emitting region for both models is also depicted in the right panel in black. The light blue line in the left and middle panels represents the H2O snow surface.

In the text

	Fig. 2 Synthetic H2O (blue), 12CO2 (green), and 13CO2 (purple) spectra. The top row depicts the model using the chemical network, and the bottom row depicts the model with parameterized abundances. In the right panels, blue triangles indicate H2O lines with Eup < 2500 K.
In the text

Fig. 3 Temperature as a function of radius within the 70% emitting region of the H2O 177,10−164,13 (6371 K), 113,9−100,10 (2438 K), and 83,6−70,7 (1447 K) lines (solid red, pink, and blue lines) for the model using the chemical network (left panel) and the model with parameterized abundances (right panel). The shaded regions represent the minimum and maximum temperature within the emitting region. The red, pink, and blue crosses represent the retrieved T and Req from the single-temperature slab fits in the 10–14, 13.5–17.5, and 21–24 μm regions. The white plus symbols represent the retrieved T and Req values from the three-temperature-component MCMC routine for H2O.

In the text

Fig. 4 Vertically integrated H2O column density as a function of radius. The solid lines show the total model column density, and the dashed lines show the model column density integrated up to the dust τ = 1 surface at 15 μm. The red, pink, and blue crosses represent the retrieved N and Req from the single-temperature slab fits in the 10–14, 13.5–17.5, and 21–24 μm regions. The plus symbols in the left column represent the retrieved N and Req values from the three-temperature-component MCMC routine. The left panel depicts the model using the chemical network and the right panel depicts the model with parameterized abundances.

In the text

	Fig. 5 Synthetic H2O spectra of the fiducial model using the chemical network (black) and three models with parameterized abundances (blue, pink, and red) between 23.7 and 24 μm. The flux is normalized to the 98,1−87,2 line at 23.87 μm.
In the text

Fig. 6 H2O temperature diagnostic diagram, as defined in Banzatti et al. (2025). The observations presented in Banzatti et al. (2025), Temmink et al. (2025), and Smith et al. (2025) are shown as gray plus symbols, blue crosses, and pink triangles, respectively. The parameterized models (denoted as PA) presented in this work are shown as blue diamonds, and the chemistry models are shown as red stars. In both cases, the fiducial model is indicated as the largest symbol with a thicker black outline.

In the text

	Fig. 7 H2O abundance maps of four models using the chemical network with varying accretion luminosity. The τ = 1 surface of the dust continuum at 15 μm is indicated with a dashed blue line. Red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371K), 113,9−100,10 (Eup = 2438K), and 83,6−70,7 (Eup = 1447K) lines, respectively.
In the text

	Fig. 8 Left: Input blackbody spectra for an L* = 0.8L⊙ source with an accretion luminosity of Lacc = 0, 0.002, 0.02 and 0.2L⊙ (black, yellow, pink, and blue). Right: Synthetic H2O spectra corresponding to the four models in the left panel.
In the text

Fig. 9 H2O photodissociation timescale and turbulent diffusion mixing timescale of the fiducial model using the chemical network. The 1,10 and 100 year timescale contours are indicated in white in both panels. The dark blue contours in both panels indicate the 10−5, 10−6, and 10−8 H2O abundance contours of the fiducial chemistry model (see also top left panel of Fig. 1). The light blue line represents the H2O snow surface.

In the text

	Fig. 10 Synthetic H2O spectra of the fiducial model using the chemical network (black) and the fiducial parameterized model (abundance above the snow surface of 10−4; blue) between 23 and 33, 40 and 55, 60 and 80, and 120 and 185 μm. The flux is normalized to the 98,1−87,2 line (2892 K) at 23.87 μm. Blue triangles indicate H2O lines with Eup < 2500 K.
In the text

Fig. A.1 H2O abundance maps of the fiducial model using the chemical network (first panel) and three models with parameterized abundances of 10−4, 10−5, and 10−6 beyond 1 au (second, third, and fourth panels). The τ = 1 surface of the dust continuum at 15 μm is indicated with a dashed blue line. Red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371 K), 113,9−100,10 (Eup = 2438 K), and 83,6−70,7 (Eup = 1447 K) lines, respectively.

In the text

Fig. A.2 Abundance maps of CO2 for the fiducial models with fℓ = 0.9. The left panel depicts the model using the chemical network, and the right panel depicts the model with parameterized abundances. The τ = 1 surface of the dust continuum at 15 μm is indicated with a dashed blue line in all panels. The green and purple contours represent the 70% emitting regions of the 12CO2 and 13CO2 ν2 = 1–0 Q(20) lines, respectively, and the pink contours represent the 70% emitting region of the H2O 113,9−100,10 (2438 K) line.

In the text

Fig. A.3 Temperature as a function of radius within the 70% emitting region of the 12CO2 and 13CO2 ν2 = 1–0 Q(20) lines (solid green and purple lines). The shaded regions represent the minimum and maximum temperature within the emitting region. The green and purple crosses represent the retrieved T and Req for 12CO2 and 13CO2. The left panel depicts the model using the chemical network, and the right panel depicts the model with parameterized abundances.

In the text

Fig. A.4 Vertically integrated CO2 column density as a function of radius. The solid lines show the total model column density, and the dashed lines show the model column density integrated up to the dust τ = 1 surface at 15 μm. The green and purple crosses represent the retrieved N and R for 12CO2 and 13CO2. The left panel depicts the model using the chemical network, and the right panel depicts the model with parameterized abundances.

In the text

Fig. A.5 Temperature as a function of radius within the 70% emitting region of the H2O 177,10−164,13 (6371 K), 113,9−100,10 (2438 K), and 83,6−70,7 (1447 K) lines (solid red, pink, and blue lines) for the model using the chemical network (left panel) and the model using parameterized abundances (right panel) with the excitation calculated in LTE. The red, pink, and blue crosses represent the retrieved T and Req from the single-temperature slab fits in the 10-14, 13.5-17.5, and 21-24 μm regions. The opaque crosses represent the fits done on the LTE-excitation models, and the faded crosses represent those done on the non-LTE excitation models (corresponding to the fits shown in Fig. 3).

In the text

Fig. A.6 Temperature as a function of radius within the 70% emitting region of the H2O 177,10−164,13 (6371 K), 113,9−100,10 (2438 K), and 83,6−70,7 (1447 K) lines (solid red, pink, and blue lines) for the models with fℓ = 0.999 with the full chemistry (left panel) and with parameterized abundances (right panel). The shaded regions represent the minimum and maximum temperature within the emitting region. The red, pink, and blue crosses represent the retrieved T and Req from the single-temperature slab fits in the 10-14, 13.5-17.5, and 21-24 μm regions. The white plus symbols represent the retrieved T and Req values from the three-temperature-component MCMC routine.

In the text

Fig. A.7 Vertically integrated H2O column density as a function of radius for the models with fℓ = 0.999 with the full chemistry (left panel) and with parameterized abundances (right panel). The solid lines show the total model column density, and the dashed lines show the model column density integrated up to the dust τ = 1 surface at 15 μm. The red, pink, and blue crosses represent the retrieved N and Req from the single-temperature slab fits in the 10-14, 13.5-17.5, and 21-24 μm regions. The white plus symbols represent the retrieved N and Req values from the three-temperature-component MCMC routine.

In the text

	Fig. A.8 Maps of the H2 density, UV radiation field G0, and gas temperature of the fiducial chemistry model. The solid black (and white in the second panel) contours represent the H2O abundance structure of the fiducial parameterized model (see Fig. 1).
In the text

	Fig. A.9 H2O abundance as a function of gas temperature and UV radiation field, following a simple chemical model (see Sect. A.4). The red and blue boxes represent the conditions within the inner and outer H2O reservoir.
In the text

Fig. A.10 H2O abundance, gas temperature, and H2O spectrum between 23.7 and 24 μm for a model with a gas temperature increased by a factor of 2 with respect to the fiducial chemistry model. In the first panel, red, pink, and blue contours represent the 70% emitting regions of the H2O 177,10−164,13 (Eup = 6371K), 113,9−100,10 (Eup = 2438K; also shown in black in the second panel), and 83,6−70,7 (Eup = 1447K) lines, respectively. The dust τ = 1 surface at 15 μm is indicated with a dashed blue line in the first and second panel. In the third panel, the fiducial chemistry (chem.) model is shown in black; the model with increased temperature is shown in red; and the fiducial parameterized model (denoted as PA) is shown in blue.

In the text