# Understanding the water emission in the mid- and far-IR from protoplanetary disks around T Tauri stars

S. Antonellini, I. Kamp, P. Riviere-Marichalar, R. Meijerink, P. Woitke, W.

-F. Thi, M. Spaans, G. Aresu, E. Lee

► To cite this version:

S. Antonellini, I. Kamp, P. Riviere-Marichalar, R. Meijerink, P. Woitke, et al. Understanding the water emission in the mid- and far-IR from protoplanetary disks around T Tauri stars. Astronomy and Astrophysics - A&A, 2015, 582, 10.1051/0004-6361/201525724. insu-03634052

HAL Id: insu-03634052

https://insu.hal.science/insu-03634052

Submitted on 7 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Understanding the water emission in the mid- and far-IR from protoplanetary disks around T Tauri stars*

S. Antonellini<fnref n="1" />, I. Kamp<fnref n="1" />, P. Riviere-Marichalar<fnref n="2" />, R. Meijerink<fnref n="3" />, P. Woitke<fnref n="4" />, W.-F. Thi<fnref n="5" />,<fnref n="6" />, M. Spaans<fnref n="1" />, G. Aresu<fnref n="1" />,<fnref n="7" />, and E. Lee<fnref n="4" />

Received 22 January 2015 / Accepted 19 July 2015

## ABSTRACT

**Aims.** We investigate which properties of protoplanetary disks around T Tauri stars affect the physics and chemistry in the regions where mid- and far-IR water lines originate and their respective line fluxes. We search for diagnostics for future observations.

**Methods.** With the code ProDiMo, we build a series of models exploring a large parameter space, computing rotational and ro-vibrational transitions of water in nonlocal thermodynamic equilibrium (non-LTE). We select a sample of transitions in the mid-IR regime and the fundamental ortho and para water transitions in the far-IR. We investigate the chemistry and the local physical conditions in the line emitting regions. We calculate *Spitzer* spectra for each model and compare far-IR and mid-IR lines. In addition, we use mid-IR colors to tie the water line predictions to the dust continuum.

**Results.** Parameters affecting the water line fluxes in disks by more than a factor of three are: the disk gas mass, the dust-to-gas mass ratio, the dust maximum grain size, interstellar medium (ISM) UV radiation field, the mixing parameter of Dubrulle settling, the disk flaring parameter, and the dust size distribution. The first four parameters affect the mid-IR lines much more than the far-IR lines.

**Conclusions.** A key driver behind water spectroscopy is the dust opacity, which sets the location of the water line emitting region. We identify three types of parameters, including those (1) affecting global disk opacity and opacity function (maximum dust size and dust size distribution); (2) affecting global disk opacity (dust-to-gas mass ratio, Dubrulle settling, disk gas mass); and (3) not affecting disk opacity (flaring parameter, ISM UV radiation field, fraction of PAHs). Parameters, such as dust-to-gas ratio, ISM radiation field, and dust size distribution, affect the mid-IR lines more, while the far-IR transitions are more affected by the flaring index. The gas mass greatly affects lines in both regimes. Higher spectral resolution and line sensitivities, like from the *James Webb* Space Telescope, are needed to detect a statistically relevant sample of individual water lines to distinguish further between these types of parameters.

**Key words.** protoplanetary disks – line: formation – stars: pre-main sequence

## 1. Introduction

Water is a fundamental component of protoplanetary disks. Its role is already pivotal in the disk formation phase as an important coolant in the cloud collapse (see, e.g., recent works, Zinnecker & Yorke 2007; van der Tak et al. 2013; van Dishoeck et al. 2014). Water ice mantles can facilitate the sticking between dust grains (Machida & Abe 2010). To date water has even been detected in several exoplanets (e.g., Brogi et al. 2014; Fraine et al. 2014; Crouzet et al. 2014), and is common throughout the whole solar system in comets, meteorites, moons, and rocky planets (Encrenaz 2008). The observed abundance of water on Earth, Mars, and maybe Venus (combined with D/H ratio measurements) suggests a delivery process by comets (van Dishoeck et al. 2014).

Water in protoplanetary disks around T Tauri stars has been detected over a wide spectral range: near-IR (NIR), ground-based spectra (Salyk et al. 2008), mid-IR *Spitzer*, VISIR spectra

(Carr & Najita 2008; Pontoppidan et al. 2010a,b), and far-IR Herschel spectra (Hogerheijde et al. 2011; Riviere-Marichalar et al. 2012). For *Spitzer*, the detected transitions span from 5 to 36 µm, and they originate from warm water gas layers (T ≈ 100–2000 K) close to the planet-forming region. At the low resolution of *Spitzer*, these lines are blended in complexes with very close upper level energies (*E*<sub>up</sub>). It is possible to use different *E*<sub>up</sub> of the transitions to probe the radial distribution of water in disks (Pontoppidan et al. 2010b; Zhang et al. 2013). *Spitzer* observations found a correlation between water and HCN in T Tauri stars. The presence of these molecular features can be explained with low continuum opacity (possibly due to settling; Najita et al. 2013). Riviere-Marichalar et al. (2012) reported Herschel/PACS observations of o-H<sub>2</sub>O 63.32 µm (8<sub>18</sub>−7<sub>07</sub>, *E*<sub>up</sub> = 1071 K). The Herschel/HIFI water lines have only been detected for TW Hya (Hogerheijde et al. 2011) and DG Tau (Podio et al. 2013). Bergin et al. (2010) report a tentative detection toward DM Tau. The reason for this low detection rate can be due to water first freezing on dust, then settling and possibly ice transportation to the midplane beyond the

<footnote n="1">¹ Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands
e-mail: antonellini@astro.rug.nl</footnote>

<footnote n="2">² Centro de Astrobiología (INTA-CSIC) – Depto. Astrofísica, POB 78, ESAC Campus, 28691 Villanueva de la Cañada, Spain</footnote>

<footnote n="3">³ Leiden Observatory, Leiden University, PO Box, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4">⁴ SUPA, School of Physics and Astronomy, University of St. Andrews, St. Andrews KY16 9SS, UK</footnote>

<footnote n="5">⁵ Université Joseph Fourier Grenoble-1, CNRS-INSU, Institut de Planétologie et d'Astrophysique (IPAG) UMR 5274, 38041 Grenoble, France</footnote>

<footnote n="6">⁶ Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="7">⁷ INAF, Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius, Italy</footnote>

<footnote>* Appendices are available in electronic form at http://www.aanda.org</footnote>
snow line (Meijerink et al. 2009; Hogerheijde et al. 2011, 2012), or because of photodissociation in the outer disk (Dominik et al. 2005). Two modeling studies for TW Hya and DG Tau performed with ProDiMo show different gas temperatures in these regions (≈30 K for TW Hya; Kamp et al. 2013; ≈50–600 K for DG Tau; Podio et al. 2013).

Next to these observational studies, disk modeling has been carried out using various degrees of complexity and self-consistency in the numerical codes. Meijerink et al. (2008) modeled far-IR water transitions with an X-ray thermochemical code using a typical T Tauri disk model. They show that the far-IR lines probe the physical conditions in the outer disk (10–100 AU) and the level of turbulence. These lines are typically in non-LTE because of their large critical densities ($n_{cr} > 10^{8–12} \text{ cm}^{-3}$ Meijerink et al. 2009). Using a simplified water abundance prescription, Meijerink et al. (2009) find that dust depletion (due to settling, dust growth, lower dust-to-gas ratio) is required to produce the observed line-to-continuum ratio in the far-IR line emitting regions. Glassgold et al. (2009), using an X-ray irradiated chemical slab model at 1 AU, reproduced molecular column densities of water consistent with the observations of IR lines toward T Tauri disks when H₂ is efficient formation on dust grains is considered. A DM Tau study performed with a 1+1D disk model and very simple chemical network (Dominik et al. 2005) found an extended cold water reservoir (up to the outer disk, 800 AU). This reservoir is unaffected by the external far-UV (FUV) radiation field because of a balance between photodesorption and photodissociation. One of the main formation channels for water is through ion molecule chemistry. The particular disk model for DM Tau predicted the sub-mm water lines to be in absorption and undetectable by the Herschel Space Telescope because of beam dilution.

Woitke et al. (2009b) used a radiation thermochemical disk model for a typical disk around a Herbig star to identify the presence of three water reservoirs: an inner one in LTE ($\log_e(H_2O) \sim 10^{-4}$; $\epsilon(H_2O) = n_{H_2O}/n_{tot}$), extended up to 10 AU dominated by equilibrium chemistry and endothermic reactions; a cold outer belt ($R \sim 20–150$ AU; $z/r \le 0.05$) with desorption and adsorption in tight equilibrium, where the water vapour is mainly photodissociated by the stellar FUV photons; and an intermediate irradiated region ($R \sim 1–30$ AU; $z/r \le 0.1–0.3$), where water is mainly produced by neutral-neutral channels and destroyed by photodissociation. Far-IR rotational lines with higher $E_{up}$ are produced from surface layers ($z/r \sim 0.1–0.3$). Woitke et al. also found that a computation of the level populations using vertical escape probability underestimates the population compared to a detailed Monte Carlo computation, affecting the final line fluxes by 2–45%.

Meijerink et al. (2009) found that the role of X-ray heating (10–40% efficiency) for neutral-neutral water formation is more important than FUV heating. If X-ray heating dominates locally, then far-IR lines are stronger. Simultaneously, also ion-molecule reactions become more important (Meijerink et al. 2012).

Recent modeling by Heinzeller et al. (2011) performed with a thermochemical code found that mixing and settling are important in limiting the water formation if the high H₂ formation rate on dust grains from Cazaux & Tielens (2002) is used. Including these effects, they find a better agreement with the observed column densities. Non-LTE is required to reproduce the observed mid-IR line fluxes.

Du & Bergin (2014) use radiation thermochemical disk models to study the impact of the stellar radiation field and disk properties on the water column densities/masses in a typical T Tauri disk with a focus on the inner warm reservoir (≤5 AU). At very

low dust opacities (dust-to-gas mass ratio <10⁻⁴), photodissociation can suppress the gas-phase water reservoir inside 1.3 AU. Large dust grains with a smaller scale height than the small grain population (mimicking settling) produces a similar effect.

What is still lacking after the history of modeling described above is a detailed systematic study of how the disk and stellar parameters affect the key water emission lines in the mid- and far-IR. We want to understand what conclusions on specific pure disk properties can be drawn, given the existing *Spitzer* and *Herschel* observations, and also in the light of the new capabilities with JWST/MIRI. This kind of study requires a combination of radiation thermochemical disk models with non-LTE treatment of water level populations, thus taking the dust for the underlying physical disk structure into account, as well as the chemical structure for the subsequent computation of emission lines (fraction of water column above the dust continuum). This is the aim of the work presented here. In a next paper (Antonellini et al., in prep.), we will address the effect of a different central star on the disk mid-IR water spectroscopy, trying to find the reason for the different behavior in observations toward low-mass and high-mass pre-main-sequence stars.

In Sect. 2, we describe our code, the line selection, the standard disk, and the series of models, also motivating the explored parameter space. In Sect. 3, we present the results of the standard model, its chemistry, and water spectroscopy. In Sect. 4, we present the results of our parameter study, and finally, in Sects. 5 and 6, we present the discussion and conclusions.

## 2. Modeling procedure

Our code ProDiMo (Protoplanetary Disk Model) is a radiation thermochemical code, which self-consistently computes the physical and chemical structure of disks. The code includes 2D radiative transfer, and gas chemistry, including freeze-out and desorption processes; it can perform hydrostatic or parametric description of an axisymmetric disk structure (Woitke et al. 2009a). X-ray chemistry has been implemented by Aresu et al. (2011). The code computes the energy balance of the gas, including photoelectric heating, C photoionization, H₂ photodissociation and formation heating, collisional de-excitation of H₂, cosmic ray heating, viscous heating, cooling via Ly-α, [OI] 6300 Å, and more than 10 000 atomic/molecular emission lines. Thermalization between dust and gas and free-free processes complete the energy balance.

The radiation field comprises the central star and the ISM (InterStellar Medium) UV field. The multiwavelength radiative transfer is performed using accelerated Λ iteration, and ray tracing (generally 100 rays). The statistical equilibrium of atoms/molecules, however, is computed considering a two-directional escape probability, assuming that the line source function varies slowly in local environments, where the optical depths grows quickly, and photon escape is mainly in the vertical direction. For each species treated in the spectroscopy, the code ProDiMo computes the statistical equilibrium based on a large set of molecular and collisional data. Optionally, a detailed line radiative transfer can be performed for a selected subsample of transitions.

The code has been updated with several features: a detailed computation of the photoionization and photodissociation cross-sections from van Dishoeck et al. (2006), an extension of the collisional partners for the water statistical equilibrium (Kamp et al. 2010), a soft-edge truncation of the outer disk (Sect. 2.1.3), polycyclic aromatic hydrocarbons (PAH) ionization balance and heating/cooling, UV fluorescence pumping, a
parametric description of settling, pumping by OH photodissociation, H₂ pumping by formation on dust grains, and chemical heating (Woitke et al. 2011). A prescription for the settling following the formalism of Dubrulle et al. (1995) has been implemented and is described in Sect. 2.1.2.

The chemical modeling has been performed using the UMIST 2006 database (Woodall et al. 2007). Even though the UMIST 2012 database (McElroy et al. 2013) is now available, it was not implemented when our work started. If three-body reactions are included, however, the results for water are comparable with the older network (Kamp et al., in prep.). A test performed with the UMIST 2012 database and the standard model shows that the disk gaseous water distribution is unchanged (except a tiny abundance variation, which is very difficult to quantify in the outer disk), and the line fluxes of the far-IR water lines are affected by less than 25% and those of the mid-IR lines by about 1%. Our chemical network contains a complete series of gas-phase reactions for 94 atomic, ionic, and molecular species, involving water and related precursors. Surface reactions for water are not included in the network, as they are not relevant for the emitting regions of the transitions that we consider (Kamp et al. 2013). The chemistry is calculated assuming steady state. The chemical network and heating/cooling balance are solved iteratively. Gas heating via exothermic reactions is approximated via the chemical heating efficiency (Woitke et al. 2011) using the reaction enthalpies at the reference temperature of 0 K from NIST<fnref n="1" /> and Millar et al. (1997). Dependence of the reaction thermodynamics on the real temperature is neglected.

We calculate the detailed water statistical equilibrium using a subset of the first 500 levels, including high-J rotational and ro-vibrational levels. The ro-vibrational levels are computed from the work of Tennyson et al. (2001) and Barber et al. (2006); the complete set of levels (824 in total) include 411 o-H₂O and 413 p-H₂O ro-vibrational levels, in which rotational levels up to J = 25 (and projections till K<sub>a,c</sub> = 18) and vibrational levels up to v₁ = 2, v₂ = 2, v₃ = 1 are considered. Our water energy levels, lines, and collisional data (Green et al. 1993; Faure & Josselin 2008; Rothman et al. 2013) are the same as adopted in Kamp et al. (2013).

The analysis of each model is based on the averaged physical and chemical quantities computed in the line emitting region that produces 50% of the radially and vertically integrated flux (for details see Appendix B).

## 2.1. Standard model and description of the parameter study

The standard model is a parametrized “monolithic” disk (properties defined in Table 1 are the same throughout the whole disk structure) surrounding a typical T Tauri K-type star with the parameters reported in Table 1. It represents a typical disk around a T Tauri star without viscous heating ($\alpha_{vis} = 0.0$), with X-ray irradiation from the central star, a UV excess, a 20% chemical heating efficiency due to exothermic reactions, a settling prescription for the dust according to Dubrulle et al. (1995) (Sect. 2.1.2), and a tapering off outer edge (Sect. 2.1.3).

We explore individual parameters around this standard model, considering a series of models in which we vary a single parameter at a time. A list of the considered parameters in the series and adopted values are reported in Table 2. For each parameter, the range spans from a factor of a few up to orders of magnitude around the initial standard value. We choose to model single parameters at a time to disentangle and understand

<table><caption>**Table 1.** Overview of the fixed model parameters for the standard T Tauri model.</caption><caption>Central star and radiation field parameters</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Photospheric temperature</td><td><i>T</i><sub>eff</sub> [K]</td><td>4400</td></tr><tr><td>Stellar mass</td><td><i>M</i><sub>*</sub> [M<sub>☉</sub>]</td><td>0.8</td></tr><tr><td>Stellar luminosity</td><td><i>L</i><sub>*</sub> [L<sub>☉</sub>]</td><td>0.7</td></tr><tr><td>FUV excess</td><td><i>L</i><sub>UV</sub>/<i>L</i><sub>*</sub></td><td>0.01</td></tr><tr><td>UV powerlaw exponent</td><td><i>p</i><sub>UV</sub></td><td>0.2</td></tr><tr><td>X-ray luminosity</td><td><i>L</i><sub>X</sub> [erg/s]</td><td>10<sup>30</sup></td></tr><tr><td>X-ray minimum energy</td><td><i>E</i><sub>min,X</sub> [keV]</td><td>0.1</td></tr><tr><td>X-ray temperature</td><td><i>T</i><sub>X</sub> [K]</td><td>10<sup>7</sup></td></tr></tbody></table>

<table><caption>Disk parameters of the standard model fixed in the series</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Radial × vertical grid points</td><td><i>N</i><sub>xx</sub> × <i>N</i><sub>zz</sub></td><td>70 × 70</td></tr><tr><td>Outer radius</td><td><i>R</i><sub>out</sub> [AU]</td><td>300</td></tr><tr><td>Inner radius</td><td><i>R</i><sub>in</sub> [AU]</td><td>0.1</td></tr><tr><td>Surface density power law index</td><td><i>ε</i></td><td>1.0</td></tr><tr><td>Minimum dust size</td><td><i>a</i><sub>min</sub> [µm]</td><td>0.05</td></tr><tr><td>Reference radius</td><td><i>R</i><sub>0</sub> [AU]</td><td>0.1</td></tr><tr><td>Scale height at reference radius</td><td><i>H</i><sub>0</sub> [AU]</td><td>3.5 × 10<sup>-3</sup></td></tr><tr><td>Tapering-off radius</td><td><i>R</i><sub>taper</sub> [AU]</td><td>200</td></tr><tr><td>Chemical heating efficiency</td><td>—</td><td>0.2</td></tr><tr><td>Settling description</td><td>—</td><td>Dubrulle</td></tr><tr><td>Cosmic ray ionization rate</td><td><i>ζ</i><sub>CRs</sub> [s<sup>-1</sup>]</td><td>1.7 × 10<sup>-17</sup></td></tr><tr><td>Distance</td><td><i>d</i> [pc]</td><td>140</td></tr><tr><td>Turbulence viscosity coefficient</td><td><i>α</i><sub>vis</sub></td><td>0</td></tr><tr><td>Disk inclination</td><td>[deg]</td><td>30</td></tr></tbody></table>

individual effects on both water spectroscopy and continuum fluxes/spectral index indicators in the same wavelength regime. Based on our results from model grids, such as DENT (Woitke et al. 2011; Kamp et al. 2011) and the beta2 grid (Woitke 2012), the interplay of parameters is often well represented by combining the effects of the individual parameters.

The series deliberately extends beyond real observed disks to facilitate the detection of correlations and to disentangle the effects of the various individual parameters on the water lines. In the following subsections, we provide a short description of the different parameters investigated in this modeling.

Several other parameters have been studied as well but we do not discuss them because they do not significantly affect the mid-IR lines, such as the cosmic ray ionization rate (less than 1% variation), the chemical heating efficiency (factor less than 3), the disk outer radius (4% variation), the radioactive decay rate of unstable nuclei and related heating (less than 1%), the tapering-off radius (about 10%), and the surface density power law index (about 10%). Some of these parameters particularly affect the optically thick regions of the disk, and will be discussed in a future paper on ice reservoirs.

### 2.1.1. Dust-to-gas mass ratio

In the intestellar medium, dust represents approximately 1% of the total mass (Flower & Pineau des Forets 1995), while this value can be different in environments in which dust grows, settles and the gas photoevaporates (as in protoplanetary disks). For example, an extreme ratio can occur in the case of a debris disk, in which the disk is essentially composed of dust, probably of secondary origin (recent review by Wyatt 2008). A common hypothesis is that the dust-to-gas mass ratio is related to the age of the object, and expected to increase with age. However, this

<footnote n="1">¹ http://www.nist.gov/</footnote>
<table><caption>**Table 2.** Disk parameters varied in the series.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Series values</th></tr></thead><tbody><tr><td>Flaring power law index</td><td>β</td><td>0.8, 0.85, 0.9, 0.95, 1.0, 1.05, 1.1, <b>1.13</b>, 1.15, 1.2, 1.25</td></tr><tr><td>Dust-to-gas mass ratio</td><td>d/g</td><td>0.001, <b>0.01</b>, 0.1, 1, 10, 100</td></tr><tr><td>Gas mass</td><td>M<sub>gas</sub> [M<sub>☉</sub>]</td><td>10<sup>-5</sup>, 10<sup>-4</sup>, 0.001, <b>0.01</b>, 0.05, 0.1</td></tr><tr><td>Mixing parameter</td><td>α<sub>set</sub></td><td>10<sup>-5</sup>, 10<sup>-4</sup>, 10<sup>-3</sup>, 10<sup>-2</sup>, <b>0.05</b>, 0.1</td></tr><tr><td>ISM radiation field</td><td>χ<sub>ISM</sub> [G<sub>0</sub><sup>*</sup>]</td><td>0.5, <b>1</b>, 10, 100, 1000, 10<sup>4</sup>, 10<sup>5</sup>, 10<sup>6</sup></td></tr><tr><td>Fraction of PAH<sup>**</sup></td><td>f<sub>PAH</sub></td><td>0.0001, 0.001, <b>0.01</b>, 0.1, 1.0</td></tr><tr><td>Maximum dust size</td><td>a<sub>max</sub> [μm]</td><td>250, 400, 500, 700, <b>1000</b>, 2000, 5000, 10<sup>4</sup>, 10<sup>5</sup></td></tr><tr><td>Dust size power law index</td><td>a<sub>pow</sub></td><td>2.0, 2.5, 3.0, <b>3.5</b>, 4.0, 4.5</td></tr></tbody></table>

**Notes.** Bold numbers refer to standard model values; (*) Draine units $G_0 = 1.6 \times 10^{-3}$ erg cm$^{-2}$ s$^{-1}$ (Draine 1978); (**) relative to ISM abundance (~3 × 10$^{-7}$ per H nucleus; Tielens 2008).

picture might be oversimplified, since there are detected objects in millimeter observations with low dust content, but still high accretion rates, for example, toward star formation environments like Cep OB2 (Sicilia-Aguilar et al. 2011). This suggests that the dust content of the disk alone cannot be considered a diagnostic for the disk age.

In our series of models, we explore a range of dust-to-gas mass ratio spanning from $10^{-3}$ up to 100.0 to reproduce the situation of dust poor and dust-rich protoplanetary disks. The highest dust-to-gas mass ratio represents the most evolved systems.

### 2.1.2. Turbulent viscosity and settling

A description of the dust settling has been implemented in ProDiMo following the approach of Dubrulle et al. (1995). The description of the adopted settling is based on the balance between gravity and vertical mechanical gas mixing effects in disks on dust grains. Since the effect is grain size dependent, it produces a different corresponding scale height for each dust grain size. Dubrulle et al. (1995) describe settling as a function of the grain size, density, relative dust and gas velocity, and the mixing parameter ($\alpha_{set}$). The prescription of the gas vertical mixing is conceptually the same as the turbulent motion responsible for the angular momentum redistribution through the disk and accretion, i.e.,

$$v = \alpha_{vis} c_s H. \qquad (1)$$

Here, $v$ is the kinematical viscosity of the gas, $c_s$ is the sound speed, and $H$ the gas scale height (Shakura & Sunyaev 1973). The two $\alpha$'s are physically the same quantities, but here we model passive disks (without viscous heating and accretion, and so $\alpha_{vis} = 0$) and use $\alpha_{set}$ to describe the amount of settling in the disk. The relation between settling, viscous accretion, and UV excess is still not fully understood and additional factors, such as stellar activity and dead zones, possibly contribute as well.

We consider a range of $\alpha_{set}$ that spans from $0.1$ to $10^{-5}$. The dust settling changes the corresponding dust scale height with respect to the gas scale height. As a consequence, settling produces a local variation of the dust-to-gas mass ratio. In particular, a low value of $\alpha_{set}$ means that the dust is less homogeneously mixed with the gas and hence strongly settled.

### 2.1.3. Structure of the disk

The geometrical shape of the disk is determined by the flaring index $\beta$ and the tapering-off radius $R_{taper}$. The scale height $H$

of the disk is parametrized as a power law function of the radius $r$

$$H(r) = H_0 \left( \frac{r}{R_0} \right)^\beta. \qquad (2)$$

The surface density prescription follows that often used in observational studies (de Gregorio-Monsalvo et al. 2013; Guilloteau et al. 2013), i.e.,

$$\Sigma(r) = \int_{0}^{z_{\max}} \rho(r, z) dz \propto r^{-\epsilon} e^{-\left(\frac{r}{R_{\text{taper}}}\right)^{2-\gamma}}. \qquad (3)$$

It includes an exponential radial cut-off (beyond the tapering-off radius $R_{taper}$), which affects the flaring structure in the outer disk. The parameter $\epsilon$ is the power law index of the surface density; $\gamma$ is an exponent that regulates the radial cut-off, $\gamma = \min(2, \epsilon)$; and $\rho(r, z)$ is the disk density, integrated from the midplane up to the surface. The vertical mass distribution of a disk follow hydrostatics, and the density decreases exponentially with $z$, accordingly to the radial scale height ($H$),

$$\rho(r, z) = \frac{\Sigma(r)}{\sqrt{2\pi H}} \exp\left(-\frac{z^2}{2H^2}\right). \qquad (4)$$

Observations of protoplanetary disks are often found to be consistent with the previous theoretical description, and have flared structures, with a scale height increasing with the radius (Kenyon & Hartmann 1987). HST scattered light images of disks in Taurus (Burrows et al. 1996; Stapelfeldt et al. 1998; Padgett et al. 1999) and in Orion (Smith et al. 2005) show the flared silhouette directly. Flaring is important for the thermal, ionization, and chemical structure of the disk. It is a relevant aspect of disk evolution, since the amount of irradiation of the outer regions affects photoevaporation (see recent review by Williams & Cieza 2011). The tapered disk profile can for example explain the differences in size of CO rotational emission and continuum images (e.g., Piétu et al. 2005; Isella et al. 2007; de Gregorio-Monsalvo et al. 2013). Viscous accretion models indeed predict this kind of exponential outer cut-off (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). In our series of models we explore very flat ($\beta < 1.0$) to very flared disks (up to $\beta = 1.25$) to study a wide dynamic range.

### 2.1.4. Dust grain sizes and size distributions

Dust grows in protoplanetary disks beyond sizes found in the ISM. This can be inferred from the SED slope in the mm range ($F \propto v^{\alpha_{mm}}$). In particular, the slope in the range 0.5–1 mm is
much shallower in protoplanetary disks ($\alpha_{mm} \approx 2-3$; Beckwith & Sargent 1991; Mannings & Emerson 1994; Andrews & Williams 2005, 2007) than in the diffuse ISM ($\alpha_{mm} \sim 4$; Boulanger et al. 1996). There is evidence that the growth continues during the disk lifetime, since a change in the slope has been detected for class I (median $\alpha_{mm} = 2.5$) and class II objects (median $\alpha_{mm} = 1.8$; Andrews & Williams 2005). These mm observations, however, only provide a lower limit to the dust size, since the free-free emission from an ionized stellar wind can contaminate the emission for $\lambda \gtrsim 1$ cm (Natta et al. 2004). Assuming that the 3.5 cm emission detected from TW Hya (Wilner et al. 2005) is entirely from thermal emission of dust, a bimodal dust model (Draine 2003) suggests the presence of grains up to ~10 cm, extended out at least to tens of AU to match the observations. However, growth of dust grains should proceed faster in the inner disk because of their dependence on orbital timescale and density. Hence the size distribution can also be a function of radius (e.g., Dullemond & Dominik 2005; Isella et al. 2006; Blum & Wurm 2008).

During the growth of dust particles, competitive processes of depletion start to happen, like settling, differential mixing, and fragmentation, which are all size dependent. Fragmentation can explain the presence of small grains in protoplanetary disks ($a < 100$ µm) after $10^4$ yrs (Dullemond & Dominik 2005; Zsom et al. 2011). Small grains experience a complex interplay between collisions and coagulation (Dullemond & Dominik 2008). Large dust bodies (in particular, around a meter in size) are also affected by radial drift, which can enhance local selective dust depletion (e.g., Brauer et al. 2008; Ormel et al. 2011). In the interstellar medium and in protoplanetary disks, the dust size distribution is often described as a power law in the form $n(a) \sim a^{-d_{pow}}$. The previous described effects can lead to a change in the power law distribution over the lifetime of protoplanetary disks with clear deviations from the MRN distribution (Mathis et al. 1977; Testi et al. 2014; Birnstiel et al. 2011). From the observation of the continuum in the mm, the power law index $\alpha_{mm}$ (see previous section) can be used to estimate the opacity function ($\kappa \propto \nu^\beta$), where $\beta$ is the power law index for the opacity function. From this, the dust size power law distribution can be inferred. $\beta$ is smaller than 1 for $a_{max} \gtrsim 0.5$ mm in a range of power law indices $a_{pow}$ from 2.5 to 3.5 (Natta et al. 2004; Draine 2006). Indexes that are smaller than 3.5 are to be expected as a consequence of the dust growth, while a value above 3.5 is indicative of a regime of fragmentation (Weidenschilling 1997). Many observations are consistent with indexes between 2.5 and 3.3, e.g., objects like TW Hya, CQ Tau, and HD 142666 (Natta et al. 2004). The dust size power law index is degenerate with the maximum dust size, however, values $\ge$4.0 do not fit the observations (Natta et al. 2004). To probe slightly beyond the observed range, we extended the upper limit to 4.5. The smaller the power law index, the higher the fraction of the total dust mass in big grains, thus lowering the UV opacity of the dust.

### 2.1.5. The ISM UV radiation field

Star formation mostly happens in large clusters with hundreds to thousands of stars (McKee & Williams 1997; Lada & Lada 2003) with a relevant population of O and B-type stars. These are intense sources of UV radiation, which generate local UV radiation fields that are several orders of magnitude stronger than the local UV field in the solar neighborhood. This strong radiation bath can erode the external surface regions of the disk. Examples of these objects are cocoons embedded in an ionized sheath, like those observed in Orion by the Hubble Space Telescope

(Bally et al. 2000). In these objects, the FUV photons produce a thick photon dominated region (PDR), whose exposed gas layers reach $T \ge 1000$ K.

The ISM UV radiation field has been extensively studied by Draine (1978). He defined the solar neighborhood integrated flux $G_0 = 1.6 \times 10^{-3}$ erg s$^{-1}$ cm$^{-2}$. To probe the effects of the external radiation field, we create models with $\chi_{ISM}$ spanning from less than the local solar neighborhood flux (0.5 $G_0$) up to $10^6$ $G_0$. The high limit reproduces extreme cases like the Orion Bar (Marconi et al. 1998) or active galaxies (van der Werf et al. 2010), which should also host protoplanetary disks as part of the star formation process.

### 2.1.6. PAHs in disks

From a study of the interstellar medium opacity function (Draine & Lee 1984), dust is assumed to be a mixture of silicates and carbon compounds; however, PAHs contribute as well to the total opacity function. PAHs are of major importance for the gas physics of the disks, since they represent an important source of heating in the irradiated layers (Kamp & Dullemond 2004). Usually, for protoplanetary disks, the fraction of PAHs ($f_{PAHs}$) is assumed to be 0.01 with respect to the ISM abundance (~3 $\times 10^{-7}$ per H nucleus; Tielens 2008). However, there are no clear observational constraints because for T Tauri stars the PAH features are absent with relatively few exceptions (Gürtler et al. 1999). A more detailed study concerning embedded YSOs suggests that the carriers of the PAHs are depleted by at least an order of magnitude with respect to the ISM (Geers et al. 2009).

The PAH abundance can affect lines produced in the upper layers of the disk or in optically thin regions. The model series assumes a range of $f_{PAH}$ from $10^{-4}$ to 1.0. In the current implementation, $f_{PAH}$ only affects the heating of the gas and is not treated as an additional opacity source.

### 2.1.7. Gas masses

Gas generally dominates the total mass of protoplanetary disks, at least for the youngest objects. However, its detection is more challenging than the dust, since the emission is in discrete lines and bands. Their detection often requires high spectral resolution. Emission lines can only be detected if the gas excitation temperature ($T_{ex}$) is larger than the dust temperature ($T_{dust}$).

The mass of the disk can be inferred from line observations in the (sub)mm, using tracers as CO and CO isotopologues ($^{13}$CO, $^{18}$O) (e.g., Zuckerman et al. 1995; Dent et al. 2005; Miotello et al. 2014). Kamp et al. (2011) propose a method for the gas mass estimation based on the [OI] 63 µm and CO $J = 2-1$ transitions for low-mass disks ($M_{gas} < 10^{-3} M_{\odot}$). Alternatively, sub-mm continuum observations can be used to derive a disk mass from the dust mass, making assumptions on the dust-to-gas mass ratio (Beckwith & Sargent 1991; Andrews & Williams 2005, 2007). Continuum observations toward Taurus-Auriga and $\rho$ Oph are consistent with disk masses in the range $10^{-4}-0.1 M_{\odot}$. To analyze the situation of the most extreme disks observed toward several SFRs, we investigate a range of $M_{gas}$ from $10^{-5}$ to $0.1 M_{\odot}$.

## 3. Standard model results

Our standard disk model (Fig. 1) shows three main water gas reservoirs (A, B, C) and an ice reservoir (I). A similar structure has been found previously for a generic Herbig disk (Woitke et al. 2009b). The classification of the different reservoirs
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Water reservoirs in the standard disk; the white line represents the contour of $T_{\text{dust}}$ = 100 K, the black and red dashed curved lines represent the contours of optical extinction for $A_V$ = 1 and $A_V$ = 10, respectively. The water ice reservoir is outlined by the white contour and the dashed line corresponds to $A_V$ = 10. Vertical red dashed lines and the orange rectangle delimit the water reservoirs.</figcaption></figure>

adopted here is, however, different. Because of the overall lower disk temperature compared to Herbig disks, the ice reservoir is more extended in the case of the T Tauri disks. In general, models of circumstellar disks find similar structures for the water vapour distribution, but the size of the reservoirs change according to the detailed physical conditions (van Dishoeck et al. 2014). The radial extent of each reservoir is defined here using the density and the emitting regions of the water lines in the standard model. The ice reservoir is shaped by the disk opacity. At $A_V \ge 10$, the dust temperatures and photodesorption rates are so low that H$_2$O freezes onto dust grains. The exact position of the snow line is a function of the gas temperatures and pressure (Lecar et al. 2006).

The inner zone of the disk (A) is a warm reservoir with high density (≥10⁹ cm⁻³), which produces the *Spitzer* lines. Disk (A) zone is that with the largest water column density in the disk, and is extended from the inner wall up to the snow line. The lowest abundance gas reservoir is the outer reservoir (C). In this region, the density is below 100 cm⁻³, and the optical depth is low. Hence, transitions from this cold outer region of the disk probe deeper into the disk; examples are the Herschel/HIFI lines as 538.29 µm and 269.27 µm. Reservoir C is extended from the outer end of the snow line (45 AU) up to the outer disk radius. The intermediate region above the ice is the reservoir B ($z/r \sim 0.1–0.25$); it is radially extended between the reservoirs A and C. In this region, physical conditions (densities of H$_2$, H, e⁻, temperatures of dust and gas) show large changes with radius. Many of the Herschel/PACS and Spitzer Space Telescope high excitation lines (200 < $E_{\text{up}}$ < 1000 K, not analyzed in this work) arise in this disk zone. Below this region, water is predominantly in the form of ice (I). The icy reservoir in Fig. 1 is here defined by a combination of the $T_{\text{dust}} \approx 100$ K line and the

location of the continuum opacity<fnref n="2" /> $A_V$ = 10. In the ice reservoir, the water abundances of the gas phase are depleted by 5 orders of magnitude with respect to the upper reservoir B. The ice reservoir can be indirectly probed through ice features arising from dust covered by mantles of frozen water. The 3 µm NIR ice feature has been detected by Terada et al. (2007), while McClure et al. (2012) reported the tentative detection of the 63 µm far-IR feature. Gas and dust are thermally coupled in the regions of the line emission of reservoir A and C (Fig. C.2). In reservoir B, dust and gas temperatures are partially decoupled.

In general, dust settling produces a lower local dust content in the surface layers of all models (Fig. C.2) and steepens the local dust power law distribution with respect to a completely unsettled case. The effect of settling is particularly relevant in the outer disk, where the mixing is not strong enough to stir up the dust grains. This makes the average grain size $\langle a \rangle$ close to the midplane larger with respect to an unsettled model. In the line emitting region of the 12.407 µm water line (Fig. 2) for the standard model the average dust size is 0.45 µm, while in the region of the 538.29 µm this size drops to 0.36 µm. The standard model has a high mixing parameter, and is thus close to an unsettled disk model. In the water line emitting regions of all reservoirs ($z/r \approx 0.15–0.2$), above the $A_V > 1$ layer, chemical heating and H$_2$ formation heating are the main heating processes (Fig. C.2). The cooling is dominated by H$_2$O rotational and ro-vibrational emission and CO rotational and ro-vibrational emission. Deeper in the disk, gas-grain collisions couple $T_{\text{gas}}$ and $T_{\text{dust}}$. Additional plots on the standard model properties can be found in Appendix C (Fig. C.2).

### 3.1. Water chemistry

A complete list of the main reactions involving water in the disk model is reported in Table 3. The most important reactions in the emitting region of every transition are summarized in Table 4. They are sorted by their relevance with respect to the timescale and the number of grid points in which they dominate. The kinematic reaction constants of the gas phase are expressed in a parametrized manner with three coefficients ($\alpha, \beta, \gamma$). The calculation of the various types of reaction rates, using these coefficients, is described in Woodall et al. (2007) and Woitke et al. (2009a).

The emitting regions of the *Spitzer* transitions (except the 33.510 µm line) are co-spatial (Fig. 2). Hence the main processes and chemical rates are very similar. For these transitions, the main formation channels are the neutral-neutral reactions of H$_2$ and OH, and due to the presence of X-rays, the dissociative recombination of H$_3$O⁺. For the 33.510 µm line, there is an additional neutral-neutral nitrogen chemistry channel via NH$_2$ and NO. The reaction rates (for a two-body reaction, $k_{ij} \cdot n_i \cdot n_j$ [s⁻¹]) in these dense and hot regions of the disk are above 10³ s⁻¹. Water in these regions is chemically mainly depleted through neutral-neutral reactions involving atomic hydrogen and through photodissociation.

In the outer disk (reservoir C), water is produced mainly through dissociative recombination of H$_3$O⁺ and radiative association of OH and H. Secondary formation channels are the neutral-neutral reactions involving NO and NH$_2$. In this reservoir, several processes deplete water vapor: photoionization, charge exchange with He⁺ and ion-neutral reaction with Si⁺ (driven by X-rays, cosmic rays and collisional ionization with

<footnote n="2">² $A_V$ is the optical extinction in magnitudes; *r* and *z* are the radial and vertical coordinates: $A_V(r, z) = \min(A_V(r); A_V(z); A_V(R_{\text{out}}) - A_V(r))$.</footnote>
<table><caption>**Table 3.** Main reaction channels in order or relevance for water formation and destruction in the standard disk model.</caption><thead><tr><th>Reaction number</th><th>Reaction</th><th>α [cm<sup>3</sup> s<sup>-1</sup>] or [s<sup>-1</sup>]</th><th>β [-]</th><th>γ [K]</th><th>Type</th><th>References</th></tr></thead><tbody><tr><td>1</td><td>H<sub>2</sub> + OH → H<sub>2</sub>O + H</td><td>2.05 × 10<sup>-12</sup></td><td>1.52</td><td>1740.00</td><td>Neutral-Neutral</td><td>1</td></tr><tr><td>2</td><td>NH + OH → H<sub>2</sub>O + N</td><td>3.11 × 10<sup>-12</sup></td><td>1.20</td><td>0.00</td><td>Neutral-Neutral</td><td>2, 1</td></tr><tr><td>3</td><td>NH<sub>2</sub> + OH → H<sub>2</sub>O + NH</td><td>7.78 × 10<sup>-13</sup></td><td>1.50</td><td>-230.00</td><td>Neutral-Neutral</td><td>1</td></tr><tr><td>4</td><td>NH<sub>2</sub> + NO → N<sub>2</sub> + H<sub>2</sub>O</td><td>4.27 × 10<sup>-11</sup></td><td>-2.50</td><td>331.00</td><td>Neutral-Neutral</td><td>2, 1</td></tr><tr><td>5</td><td>OH + OH → H<sub>2</sub>O + O</td><td>1.65 × 10<sup>-12</sup></td><td>1.14</td><td>50.00</td><td>Neutral-Neutral</td><td>3</td></tr><tr><td>6</td><td>OH + H<sub>2</sub>CO → HCO + H<sub>2</sub>O</td><td>2.22 × 10<sup>-12</sup></td><td>1.42</td><td>-416.00</td><td>Neutral-Neutral</td><td>1</td></tr><tr><td>7</td><td>NH<sub>2</sub> + H<sub>3</sub>O<sup>+</sup> → H<sub>2</sub>O + NH<sub>3</sub><sup>+</sup></td><td>9.70 × 10<sup>-10</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>4</td></tr><tr><td>8</td><td>H<sub>3</sub>O<sup>+</sup> + Si → SiH<sup>+</sup> + H<sub>2</sub>O</td><td>1.80 × 10<sup>-9</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>4</td></tr><tr><td>9</td><td>H<sub>3</sub>O<sup>+</sup> + SiO → SiOH<sup>+</sup> + H<sub>2</sub>O</td><td>2.00 × 10<sup>-9</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>5</td></tr><tr><td>10</td><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → H<sub>2</sub>O + H</td><td>1.08 × 10<sup>-7</sup></td><td>-0.50</td><td>0.00</td><td>Dissociative recombination</td><td>6</td></tr><tr><td>11</td><td>H<sup>-</sup> + OH → H<sub>2</sub>O + e<sup>-</sup></td><td>1.00 × 10<sup>-10</sup></td><td>-5.22</td><td>90.00</td><td>Associative detachment</td><td>4</td></tr><tr><td>12</td><td>H + OH → H<sub>2</sub>O + hν</td><td>5.26 × 10<sup>-18</sup></td><td>-5.22</td><td>90.00</td><td>Radiative association</td><td>7</td></tr><tr><td>13</td><td>H<sub>3</sub>O<sup>+</sup> + hν → H<sub>2</sub>O + H<sup>+</sup></td><td>2.00 × 10<sup>-11</sup></td><td>0.00</td><td>2.00</td><td>Photodissociation</td><td>8</td></tr><tr><td>14</td><td>H<sub>2</sub><sup>+</sup> + OH → H<sub>2</sub>O + H</td><td>3.60 × 10<sup>-11</sup></td><td>0.00</td><td>0.00</td><td>Neutral-Neutral<sup>a</sup></td><td>1</td></tr><tr><td>15</td><td>H<sub>2</sub>O<sub>(ice)</sub> + dust → H<sub>2</sub>O + dust</td><td>—</td><td>—</td><td>4800.00</td><td>Desorption<sup>b</sup></td><td>9, 19</td></tr><tr><td>16</td><td>H<sub>2</sub>O<sub>(ice)</sub> + C.R.P. → H<sub>2</sub>O + dust</td><td>—</td><td>—</td><td>4800.00</td><td>Cosmic-rays-desorption<sup>b</sup></td><td>19</td></tr><tr><td>17</td><td>H<sub>2</sub>O<sub>(ice)</sub> + hν → H<sub>2</sub>O</td><td>—</td><td>—</td><td>4800.00</td><td>Photodesorption<sup>b</sup></td><td>9, 19</td></tr><tr><td>18</td><td>PAH<sup>-</sup> + H<sub>3</sub>O<sup>+</sup> → PAH + H<sub>2</sub>O + H</td><td>1.70 × 10<sup>-8</sup></td><td>0.50</td><td>0.00</td><td>PAH dissociative recombination</td><td>18</td></tr><tr><td>1</td><td>H<sub>2</sub>O + H → OH + H<sub>2</sub></td><td>1.59 × 10<sup>-11</sup></td><td>1.20</td><td>9610.00</td><td>Neutral-Neutral</td><td>1</td></tr><tr><td>2</td><td>H<sub>3</sub><sup>+</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + H<sub>2</sub></td><td>5.90 × 10<sup>-9</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>10</td></tr><tr><td>3</td><td>C<sup>+</sup> + H<sub>2</sub>O → HCO<sup>+</sup> + H</td><td>9.00 × 10<sup>-10</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>4</td></tr><tr><td>4</td><td>CH<sub>5</sub><sup>+</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + CH<sub>4</sub></td><td>3.70 × 10<sup>-9</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>10</td></tr><tr><td>5</td><td>H<sub>2</sub>O + Si<sup>+</sup> → SiOH<sup>+</sup> + H</td><td>2.30 × 10<sup>-10</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>11</td></tr><tr><td>6</td><td>H<sub>2</sub>O + HCO<sup>+</sup> → CO + H<sub>3</sub>O<sup>+</sup></td><td>2.50 × 10<sup>-9</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>12</td></tr><tr><td>7</td><td>H<sub>2</sub>O + HN<sub>2</sub><sup>+</sup> → N<sub>2</sub> + H<sub>3</sub>O<sup>+</sup></td><td>2.60 × 10<sup>-9</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>13</td></tr><tr><td>8</td><td>H<sub>2</sub>O + SiH<sup>+</sup> → Si<sup>+</sup> + H<sub>3</sub>O<sup>+</sup></td><td>8.00 × 10<sup>-10</sup></td><td>0.00</td><td>0.00</td><td>Ion-Neutral</td><td>5</td></tr><tr><td>9</td><td>H<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + H</td><td>6.90 × 10<sup>-9</sup></td><td>0.00</td><td>0.00</td><td>Charge exchange</td><td>14</td></tr><tr><td>10</td><td>He<sup>+</sup> + H<sub>2</sub>O → He + H<sub>2</sub>O<sup>+</sup></td><td>6.05 × 10<sup>-11</sup></td><td>0.00</td><td>0.00</td><td>Charge exchange</td><td>15</td></tr><tr><td>11</td><td>H<sub>2</sub>O + hν → OH + H</td><td>5.90 × 10<sup>-10</sup></td><td>0.00</td><td>1.70</td><td>Photodissociation</td><td>16</td></tr><tr><td>12</td><td>H<sub>2</sub>O + hν → H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup></td><td>3.30 × 10<sup>-11</sup></td><td>0.00</td><td>3.90</td><td>Photoionization</td><td>16</td></tr><tr><td>13</td><td>H<sub>2</sub>O + dust → H<sub>2</sub>O<sub>(ice)</sub></td><td>—</td><td>—</td><td>—</td><td>Adsorption<sup>b</sup></td><td>19</td></tr><tr><td>14</td><td>Ne<sup>+</sup> + H<sub>2</sub>O → Ne + H<sub>2</sub>O<sup>+</sup></td><td>8.00 × 10<sup>-10</sup></td><td>0.00</td><td>0.00</td><td>Charge exchange</td><td>17</td></tr></tbody></table>

**References.** (1) Mallard et al. (1994); (2) McElroy et al. (2013); (3) Woodall et al. (2007); (4) Prasad & Huntress (1980); (5) Herbst et al. (1989); (6) Tielens & Hollenbach (1985); (7) Field et al. (1980); (8) Jensen et al. (2000); (9) Hollenbach & Gorti (2009); (10) Kim & Huntress (1975); (11) Fahey et al. (1981); (12) Adams et al. (1978); (13) Rakshit (1982); (14) Smith et al. (1992); (15) Mauclaire et al. (1978); (16) van Dishoeck et al. (2006); (17) Anicich (1993); (18) Flower & Pineau des Forêts (2003); (19) Aikawa et al. (1996).

**Notes.** (a) $H_2$ in a vibrationally excited state, populated by fluorescence from the electronically excited states; (b) For desorption processes, $\gamma$ is the adsorption energy.

<table><caption>**Table 4.** Main reaction channels in the emitting regions of the considered transitions in the standard model.</caption><thead><tr><th>Species & wavelength [µm]</th><th>Formation reaction number</th><th>Destruction reaction number</th></tr></thead><tbody><tr><td>o-H<sub>2</sub>O 12.396</td><td>1, 10</td><td>1, 11</td></tr><tr><td>p-H<sub>2</sub>O 12.407</td><td>1, 10</td><td>1, 11</td></tr><tr><td>o-H<sub>2</sub>O 12.445</td><td>1, 10</td><td>1, 11</td></tr><tr><td>o-H<sub>2</sub>O 12.453</td><td>1, 10</td><td>1, 11</td></tr><tr><td>p-H<sub>2</sub>O 12.832</td><td>1, 10</td><td>1, 11</td></tr><tr><td>p-H<sub>2</sub>O 12.894</td><td>1, 10</td><td>1, 11</td></tr><tr><td>o-H<sub>2</sub>O 33.510</td><td>1, 10, 4</td><td>1, 11</td></tr><tr><td>p-H<sub>2</sub>O 269.27</td><td>12, 10, 4</td><td>12, 10, 13, 11, 5</td></tr><tr><td>o-H<sub>2</sub>O 538.29</td><td>12, 10, 4</td><td>12, 10, 13, 11, 5</td></tr></tbody></table>

**Notes.** Reaction numbers in order of relevance.

electrons), water photodissociation, and ice formation. The reaction rates here are lower than in reservoir A, of the order of $10^{-5}$ s$^{-1}$ (~1 day), but the timescales are still shorter than the disk lifetime.

Previous studies found differences in terms of nitrogen chemistry for T Tauri disks. The channels mentioned above are an important formation path for water precursors (as OH) in the outer disk of TW Hya, and in a tiny region in the inner rim

where water is directly formed from NO and NH$_2$ (Kamp et al. 2013). In the Najita et al. (2011) chemistry network, atomic nitrogen is a competitor that depletes OH in the disk atmosphere, making it harder to produce water. However, this model only includes X-rays and not FUV and also lacks water rotational cooling. Water in the inner disk is formed through the activated channel involving OH and H$_2$, in agreement with our modeling. The Glassgold et al. (2009) model does not include FUV and
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Water lines from different reservoirs in the standard model; left: reservoir A, center: reservoir A-B, right: reservoir B-C. The blue dashed curved lines are contours for gas densities of $10^8$ and $10^{10}$ cm$^{-3}$, the yellow contours are for gas temperatures of 50 K, 150 K, 600 K. Vertical red dashed lines delimit the reservoirs.</figcaption></figure>

its network is limited to 125 reactions and 25 species, including only H, He, C, and O. This model agrees with our study on the water formation mechanism through neutral-neutral reaction of OH and $H_2$, and finds that reactions with $H^+$ are the main destruction channels of water in the upper layers.

In the modeling of Herbig disks, performed with a comparable chemical network, Woitke et al. (2009b) found that dissociative recombination is not important in the inner disk, contrary to the case of T Tauri disks. The outer disk water reservoir of Herbigs is then fueled by photodesorption processes from the ices, while the depletion is again due to photodissociation. The main difference from the chemistry in T Tauri disks is the lack of channels involving ions.

## 3.2. Water spectroscopy

Water has an extremely rich spectroscopy, making it necessary to select a subsample of lines for this study. We thus focus here on the list of transitions shown in Table 5, which are based on *Spitzer* and *Herschel/HIFI* line detections in disks around T Tauri stars. A detailed investigation of each single line regarding the chemistry is only performed for the reference model (parameters are reported in boldface in Table 2).

The mid-IR transitions originate from the warmest water reservoir (A, Fig. 2), consistent with their high upper level energies (Table 5). The two sub-mm lines are produced in the outer cold disk (reservoir C). We chose the mid-IR 12 µm lines because they are also observable from the ground with VISIR on the VLT (Pontoppidan et al. 2010a) and they are physically unblended. This last aspect is very important, because our modeling is based on individual line radiative transfer. The strongest detections in the *Spitzer* regime are the blends around 15 µm and 17 µm (Pontoppidan et al. 2010b; Carr & Najita 2011; Zhang et al. 2013), but the difference from the 12 µm lines is within a factor 10. We verified that the behavior of the single components of the blends in the Pontoppidan et al. (2010b) work is the same as that of the individual lines we discussed. In the following, we plot all the line fluxes as a function of the selected parameters. The focus of this study is to understand the discrepancy in detection between *Spitzer*/IRS and *Herschel/HIFI*, i.e., between the inner and outer disk, and to find a possible connection/explanation with the physical properties of the disks.

Figure 2 shows the observed and detected *Spitzer* lines considered here; ($E_{up} \ge 3000$ K). Since reservoir A is very compact compared to the disk size, these transitions largely overlap in their emitting volume, and as a consequence, the physical conditions are very similar for all the lines. As long as the radial and vertical temperature and density gradients in these line emitting regions are small, a simplified slab approach as employed, e.g., in Carr & Najita (2008) presents a reasonable approximation. The transitions below 12.5 µm are produced in a region that is partially optically thick. A vertical integration of the dust extinction (see Appendix B) confirms that the $\tau_{dust} = 1$ surface at the wavelength of the line is above the mean vertical height of the line emitting region, so that the lines are partially buried in the continuum. Lines from reservoir B are produced in an optically thinner region, about 50% of the emitting region is optically thin. More than 30% of the 269.27 µm and 70% of the 538.29 µm line emitting regions (reservoir C) are vertically extended above the $\tau_{dust} = 1$ contour at the respective wavelengths.

The far-IR transitions are very close to LTE because the average density of $H_2$ is close to the critical density (within an order of magnitude), and so we expect them to be sensitive to $T_{gas}$ and $n_{gas}$ variations. Our mid-IR lines, instead, are produced in a very dense region of the disk, and so they form under LTE conditions. Table 6 shows that $H_2$ is the main collisional partner for all considered transitions because the $H/H_2$ ratio is always below 0.1. Because of changes in the physical conditions of the different models, the previous statements are strictly valid only for our standard disk model.

## 4. Results

The previously discussed standard model is a starting point from which we build a series of models, modifying a single parameter per disk model. In the following, we discuss the results from this series.

The line fluxes have been treated as independent transitions, hence, as physically unblended. The results reported here are based on vertical escape probability. This approach agrees within a factor of two with a more detailed and computationally more expensive radiative transfer (Fig. 3). In particular, far-IR lines agree within 30%, while mid-IR lines agree within 70%. The lines from vertical escape probability show a systematic stronger value than detailed RT lines. The reason for this
<table><caption>**Table 5.** Water lines selected for this work.</caption><thead><tr><th>Species</th><th>Wavelength [µm]</th><th>Instrument</th><th>Reservoir</th><th>E<sub>up</sub> [K]</th><th>Quantum numbers v<sub>1</sub>v<sub>2</sub>v<sub>3</sub>JK<sub>a</sub>K<sub>c</sub></th></tr></thead><tbody><tr><td>p-H<sub>2</sub>O</td><td>12.407</td><td>Spitzer IRS</td><td>A</td><td>4945.0</td><td>000-000 16-15 3-2 13-14</td></tr><tr><td>o-H<sub>2</sub>O</td><td>12.396</td><td>Spitzer IRS</td><td>A</td><td>5781.0</td><td>000-000 17-16 4-3 13-14</td></tr><tr><td>o-H<sub>2</sub>O</td><td>12.445</td><td>Spitzer IRS</td><td>A</td><td>3629.0</td><td>000-000 11-10 8-5 3-6</td></tr><tr><td>o-H<sub>2</sub>O</td><td>12.453</td><td>Spitzer IRS</td><td>A</td><td>4213.0</td><td>000-000 13-12 7-4 6-9</td></tr><tr><td>p-H<sub>2</sub>O</td><td>12.832</td><td>Spitzer IRS</td><td>A</td><td>3243.0</td><td>000-000 10-9 8-5 2-5</td></tr><tr><td>p-H<sub>2</sub>O</td><td>12.894</td><td>Spitzer IRS</td><td>A</td><td>3310.0</td><td>000-000 12-11 5-2 7-10</td></tr><tr><td>o-H<sub>2</sub>O</td><td>33.510</td><td>Spitzer IRS</td><td>A/B</td><td>2275.0</td><td>000-000 10-9 4-3 7-6</td></tr><tr><td>p-H<sub>2</sub>O</td><td>269.27</td><td>Herschel/HIFI</td><td>B/C</td><td>53.43</td><td>000-000 1-0 1-0 1-0</td></tr><tr><td>o-H<sub>2</sub>O</td><td>538.29</td><td>Herschel/HIFI</td><td>B/C</td><td>60.96</td><td>000-000 1-1 1-0 0-1</td></tr></tbody></table>

<table><caption>**Table 6.** Overview and coefficients of the considered transitions.</caption><thead><tr><th>Species</th><th>Wavelength [µm]</th><th>A<sub>ul</sub> [s<sup>-1</sup>]</th><th>⟨n<sub>H<sub>2</sub></sub>⟩ C<sub>H<sub>2</sub>,ul</sub><sup>a</sup> [s<sup>-1</sup>]</th><th>⟨n<sub>H</sub>⟩ C<sub>H,ul</sub><sup>a</sup> [s<sup>-1</sup>]</th><th>⟨n<sub>e<sup>-</sup></sub>⟩ C<sub>e,ul</sub><sup>a</sup> [s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>p-H<sub>2</sub>O</td><td>12.407</td><td>4.217</td><td>4.489 × 10<sup>12</sup> (o and p)</td><td>2.838 × 10<sup>9</sup></td><td>8.918 × 10<sup>-3</sup></td></tr><tr><td>o-H<sub>2</sub>O</td><td>12.396</td><td>7.665</td><td>4.823 × 10<sup>12</sup> (o and p)</td><td>3.076 × 10<sup>9</sup></td><td>1.757 × 10<sup>-2</sup></td></tr><tr><td>o-H<sub>2</sub>O</td><td>12.445</td><td>2.937 × 10<sup>-1</sup></td><td>4.982 × 10<sup>11</sup> (o and p)</td><td>2.990 × 10<sup>8</sup></td><td>4.089 × 10<sup>-4</sup></td></tr><tr><td>o-H<sub>2</sub>O</td><td>12.453</td><td>1.160</td><td>1.660 × 10<sup>12</sup> (o and p)</td><td>9.963 × 10<sup>8</sup></td><td>1.431 × 10<sup>-3</sup></td></tr><tr><td>p-H<sub>2</sub>O</td><td>12.832</td><td>1.564 × 10<sup>-1</sup></td><td>1.375 × 10<sup>12</sup> (o and p)</td><td>8.253 × 10<sup>8</sup></td><td>2.399 × 10<sup>-4</sup></td></tr><tr><td>p-H<sub>2</sub>O</td><td>12.894</td><td>2.596 × 10<sup>-1</sup></td><td>5.856 × 10<sup>11</sup> (o and p)</td><td>3.154 × 10<sup>8</sup></td><td>3.956 × 10<sup>-4</sup></td></tr><tr><td>o-H<sub>2</sub>O</td><td>33.510</td><td>4.691</td><td>1.877 × 10<sup>11</sup> (o and p)</td><td>6.689 × 10<sup>7</sup></td><td>2.862 × 10<sup>-2</sup></td></tr><tr><td>p-H<sub>2</sub>O</td><td>269.27</td><td>1.852 × 10<sup>-2</sup></td><td>1.817 × 10<sup>-3</sup> (para)</td><td>1.070 × 10<sup>-6</sup></td><td>4.145 × 10<sup>-5</sup></td></tr><tr><td>o-H<sub>2</sub>O</td><td>538.29</td><td>3.477 × 10<sup>-3</sup></td><td>1.348 × 10<sup>-3</sup> (para)</td><td>5.985 × 10<sup>-7</sup></td><td>5.869 × 10<sup>-5</sup></td></tr></tbody></table>

**Notes.** (a) Collisional rates; the average densities used are computed within the emitting region of the considered line and the T<sub>gas</sub> for the collisional constants are averaged inside the emitting region (see Appendix B).

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison between the vertical escape probability flux and a detailed radiative transfer computation for the standard model. Green dashed and blue dot-dashed lines are confidence intervals of 30% and 70% from the red line, respectively.</figcaption></figure>

mismatch is twofold: the disk inclination (30 deg) causes larger column densities of gas and dust in the detailed RT calculation with respect to the vertical escape treatment. In addition, the vertical escape does not take the line self-absorption into account. The vertical escape method is a fast and powerful method to compute tens of thousands of lines for different species. Cases in which the vertical escape probability is unreliable are rare cases of absorption lines. Those are limited to very few flat models ($\beta < 0.8$), to models with $a_{\text{pow}} \ge 4.0$, or dust-to-gas mass ratio $\ge 100$, and our conclusions are not affected.

To illustrate the results, the theoretical *Spitzer* spectrum of each model is computed with a resolution of 600 (the maximum of IRS, of LH and SH modules; Houck et al. 2004). These spectra give a first qualitative result, which can be compared to observations. We also define a “mid-IR color” as the ratio of the continuum fluxes at 13.5 and 30 µm as indicator of the SED slope in the *Spitzer* wavelength range.

In the following subsections, we discuss the behavior of the selected transitions, considering the representative mid-IR and far-IR lines. We also describe in detail some average physical quantities in the line emitting regions. The computation of these quantities is described in detail in Appendix B. We compare the predicted fluxes with the sensitivity limits for *Spitzer*/IRS (5 × 10<sup>-18</sup> W/m<sup>2</sup>, 1σ, deduced from the typical upper limits; Pontoppidan et al. 2010b), Herschel/HIFI (9 × 10<sup>-19</sup> W/m<sup>2</sup>, 1 h, 5σ; de Graauw 2004), and finally JWST/MIRI (10<sup>-20</sup> W/m<sup>2</sup>, 2.7 hrs., 10 σ; Glasse et al. 2010)<sup>3</sup>. Additional plots with the 538.29 µm o-H<sub>2</sub>O line flux and line ratios 538.29 µm/12.407 µm are included in Appendix C (Figs. C.1 and C.3).

## 4.1. Dust-to-gas mass ratio

A variation of this parameter means changing the dust mass since the gas mass is fixed to 0.01 M<sub>☉</sub>.

The water total column densities in the inner disk are unaffected by the dust content (Fig. 4b). However, the *Spitzer* line at 12.4 µm is produced by a volume in which the column density is suppressed up to 14 orders of magnitude from the low to high dust-to-gas mass ratio (Fig. 4a). The emitting region of the *Spitzer* line becomes radially more extended between the low

<footnote n="3">³ http://www.stsci.edu/jwst/instruments/miri/
sensitivity/ and http://www.stsci.edu/jwst/science/
sensitivity</footnote>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

**Fig. 4.** Dust-to-gas mass ratio models. **a)** Emitting region water column density comprises 15–85% of the radially and vertically integrated flux. Vertical dotted green lines distinguish the reservoirs. **b)** Total water column densities as a function of radius. **c)** Average extension of the line emitting region that comprises 50% of the line flux radially and vertically for 12.407 µm. **d)** Average extension of the line emitting region for 538.29 µm (50% of radially and vertically integrated flux). **e)** Theoretical *Spitzer* SH/LH modules spectra (R = 600) for this model series, continuum subtracted and arbitrarily shifted. **f)** Theoretical *Spitzer* SH/LH modules spectra (R = 600) for this model series. **g)** *Spitzer* line fluxes with sensitivity limits (green dotted lines) for *Spitzer*/IRS, JWST/MIRI. **h)** Radiation field and optical extinction in the line emitting regions of 12.407 µm and 538.29 µm.
and high dust-to-gas case (Fig. 4c). This is because the total wa-
ter reservoir is unaffected, while the emitting region moves up-
ward to the disk surface with increasing dust mass (Fig. 4c); this
causes the densities in the emitting region to drop by 8 orders
of magnitude. The continuum opacity in the line emitting region
grows with increasing dust mass (Fig. 4h).

In summary, mid-IR line fluxes are driven indirectly by the
increase in dust opacity; it forces the lines to move upward in the
disk into lower density regions (Figs. 4c and h). The continuum
optical depth increase reduces the column density of emitting
water. The deviation beyond a dust-to-gas ratio of 10.0 is due to
the extremely optically thick disk, which forces the lines to be
produced from more optically thin surface layers beyond 0.5 AU.

In the outer disk, which hosts the ice reservoir, changes in the
opacity affect the photoevaporation process, increasing the total
water vapor column density (Fig. 4b). The emitting region of the
fundamental water lines gets progressively close to $A_V = 1$.
The emitting region's radial displacement with the increase of
dust mass is less pronounced than in the mid-IR case (Fig. 4c),
but the vertical displacement is significant. Also, the line flux
of the far-IR water lines is driven by an opacity effect, but to a
lesser extent compared to the mid-IR. Also, the emitting region
moves upward with increasing dust mass, however, the emitting
water column density is only reduced by two orders of magni-
tude between the most extreme models (Fig. 4a). The situation is
less extreme, because the outer disk is only marginally optically
thick (Fig. 4a).

All *Spitzer* transitions show a strong dependence on the dust-
to-gas mass ratio (Figs. 4e and f). The line-to-continuum ratio
changes by a factor 2.5 between the smallest and largest dust-to-
gas ratio, and in the dust-rich models all line features disappear
from the *Spitzer*/IRS spectra. This is in agreement with previous
results of Meijerink et al. (2009). The continuum fluxes decrease
with decreasing dust content, and the 13.5/30 µm continuum flux
ratio decreases with increasing dust content in the disk.

## 4.2. Dust settling

Settling affects the water spectroscopy much less than the dust-
to-gas mass ratio. In fact, the largest systematic difference in wa-
ter column density can be found in the outer disk, $r > 10$ AU
(Fig. A.1b). The *Spitzer* emitting region only shows a change
in the water emitting column density in the most settled model;
that is also the only one producing a slightly stronger line flux
(Fig. A.1g). The most settled model has a higher column density
in the water emitting regions of the mid-IR and far-IR lines be-
cause of the reduced continuum opacity that boosts the thermally
activated water formation channels. In the inner disk, the density
is very high, and therefore the dust is largely unsettled. In cases
of extremely low mechanical mixing only ($\alpha_{set} \le 10^{-5}$), the disk
has such a low opacity that the line emitting region moves up-
ward because of increased radial extension in the disk, and the
column density of the emitting region is higher than in an unset-
tled case (Fig. A.1b). In the emitting region of the far-IR lines,
the column density is more affected (Fig. A.1b); for $\alpha_{set} < 0.01$,
the total water column density grows about an order of mag-
nitude (Fig. A.1d). Because of the migration of the far-IR line
emitting region, dust extinction fades quickly from a marginally
optically thick regime to an optically thin regime (Fig. A.1h).
The outer disk is affected by settling only for $\alpha_{set} < 10^{-3}$; this
threshold can be found in the evolution of several properties,
such as the opacity (Fig. A.1h).

In the outer disk the photodesorption of water (Sect. 3.1)
and photodissociation processes are not strong enough to deplete

water significantly. This produces an upward migration of the
emitting region for the far-IR transition in the most settled mod-
els. The *Spitzer* line emitting region migration (Fig. A.1c) and
larger emitting water column density are negligible given nu-
merical uncertainties.

The reduction of the mixing parameter produces a mod-
erate enhancement in the mid-IR line-to-continuum ratio
(Figs. A.1e and f), along with a reduction of the continuum flux
because of the decrease of dust in the surface layers of the in-
ner disk. The 13.5/30.0 µm continuum flux ratio decreases with
increasing the mixing coefficient $\alpha_{set}$.

## 4.3. Flaring parameter $\beta$

In the parametrized modeling, the disk structure can be set
through the scale height at a reference radius $H_0$ and the flar-
ing index $\beta$. Here, we only change the flaring and fix the scale
height at the inner disk radius. A value of $\beta < 1.0$ produces a
disk collapsed vertically behind the inner radius and completely
shadowed by the vertical extent of the inner disk rim.

An increasing flaring index decreases the water column den-
sity of the emitting region by up to two orders of magnitude,
and the emitting region becomes up to a factor 1.5 more radially
extended (Fig. A.2h). As consequence of the increased radial ex-
tension, the emitting region of the *Spitzer* lines also moves up-
ward with increasing flaring index (Fig. A.2c). The line fluxes
show a minimum for $\beta = 1.0$ (Fig. A.2g) and increase for
$\beta > 1.0$ up to a factor 1.5. The mid-IR spectral changes can be
explained by a pure geometrical effect. Although the variation in
the inner disk is small (since the inner disk scale height is fixed),
it is enough to increase the heating of the surface layers of the
inner disk to the central star radiation.

The outer disk total water column density is more affected
than the inner disk, and the total amount of water (beyond
80 AU) is a factor ~2 larger in the most flared case with re-
spect to the flattest case (Fig. A.2b). With increasing flaring,
the emitting region moves upward and gets vertically more ex-
tended (Fig. A.2d) because of an increase of the line opacity. The
column density of water vapor in the emitting region becomes
smaller, but the $T_{gas}$ becomes larger. This is consistent with the
monotonic increase of the Herschel/HIFI line flux.

With increasing flaring, the disk surface is irradiated more
strongly by the central star. This increases the dust tempera-
ture through the entire disk (especially in the surface layers) and
hence the mid-IR color value decreases (Fig. A.2f).

## 4.4. Maximum grain size

The total column density of water beyond 1 AU is affected by
the maximum grain size; the models with bigger grains have a
larger water vapor column density (Fig. A.3b). The lower the
opacity, the more depleted is the ice reservoir, both vertically and
radially, thus enriching the gas phase. The water column den-
sity in the emitting region of the *Spitzer* transition also increases
(Fig. A.3a). That region moves upward and becomes more ex-
tended with increasing maximum grain size (Fig. A.3g). Since
the total dust mass remains constant with increasing $a_{max}$, the
population of small grains is reduced. As a consequence, the av-
erage continuum opacity (Fig. A.3h) in the emitting region de-
creases with the increase of maximum grain size because the
radial gas and dust temperature gradient flattens and the emit-
ting region becomes radially more extended. This is a conse-
quence of the increased FUV flux through the whole disk. The
vertical displacement is a consequence of a radially more extended emitting region because the disk is flared and the emitting region moved outward also becomes more extended vertically (Fig. A.3c). The mid-IR line fluxes increase monotonically with increasing maximum dust size because of increased emitting area and lower continuum opacity (Figs. A.3c and h).

The emitting region of the *Herschel* lines (Fig. A.3d) moves closer to the midplane of the disk for larger $a_{max}$, maintaining roughly the same radial extent and position (Fig. A.3d). This is again due to the lower dust optical depth for the models with larger grains. The water column density increase with increasing $a_{max}$ translates into larger fluxes since the lines are optically thin.

The modeled *Spitzer* spectrum (Figs. A.3e and f) shows a decreasing mid-IR dust continuum for models with larger maximum grain sizes. There is, however, no significant variation of mid-IR color (within a factor 2).

## 4.5. Dust size power law distribution

The total water column density inside 1 AU remains constant when varying the power law index of the grain size distribution, while it decreases by several orders of magnitude in the outer disk (Fig. A.4b). The water column density in the *Spitzer* line emitting region becomes lower, and the region shrinks radially by more than a factor of 3 (Figs. A.4a and c). The migration is, however, not monotonic (Fig. A.4c). This is because in the models with a power law index above 3.0, the continuum optical depth in the *Spitzer* emitting region increases strongly. Below 3.0, the optical depth stays constant ($\tau_{cont} \sim 100$). The increase in opacity causes the line emitting region to move inward, where the gas still reaches the line excitation temperatures.

In the *Herschel*/HIFI emitting region, dust and gas are thermally coupled above a threshold of $a_{pow} = 3.5$. The emitting region becomes progressively optically thick for larger $a_{pow}$, but it never reaches $A_V = 1$. The line is emitted deeper in the disk in the models with $a_{pow} = 2.0$ up to 3.0. This is due to enhanced photodesorption, which increases the water column density in the sub-mm line emitting region. The lower column densities of water vapour in the outer disk lead to lower fluxes in the sub-mm water lines (Figs. A.4a and b). The emitting region of the far-IR lines moves upward for $a_{pow} > 3.0$. In this case, the opacity is so high that water is not efficiently photodesorbed anymore and the line has to be produced higher up in the disk, where water is still in gaseous form.

The modeled *Spitzer* spectrum (Figs. A.4e and f) shows an increase of the continuum in the NIR to mid-IR range and a lower line-to-continuum ratio for the models with a larger amount of small dust (high values of $a_{pow}$). The mid-IR color variation is very high and the 13.5/30 µm continuum flux ratio decreases by more than a factor 2. The change in continuum flux is a consequence of settling (present in all the models). For larger $a_{pow}$, the disk has more smaller grains (that are not settled), and so the warm surface layers contain a higher dust mass, which contributes to the continuum. In this case, the dust consists mainly of small grains and the opacity is high. As discussed in Sect. 4.1, the line-to-continuum ratio is then very low. The flat grain size distribution model has a lower opacity (flatter and more than 4 dex weaker in the UV-optical range), and the dust and gas temperature gradients are radially and vertically steeper than in the case of a steeper grain size distribution. Then, the dust becomes colder overall, producing less mid-IR continuum.

## 4.6. ISM radiation field

The ISM FUV radiation field affects the outer regions of the disk and changes the water column density beyond 30 AU in a complex way (Fig. A.5b). In the most irradiated model, the total water column density decreases from photodissociation. In the emitting region of the *Spitzer* lines, the column density is unaffected below $\chi_{ISM} = 10^3 G_0$ (Fig. A.5a). The *Spitzer* line emitting region for $\chi_{ISM} > 10^3 G_0$ extends over the entire disk, overlapping with the sub-mm lines. This happens because of a strong change in the thermal structure of the irradiated disk, and the formation of a hot layer of gas extending into the outer disk (beyond 10 AU). The hot gas extends down to the midplane for the $\chi_{ISM} = 10^6 G_0$ model. Water is then eroded through photodissociation in the outer disk and the line is produced at smaller radii and closer to the midplane (Fig. A.5d). The increase in $T_{gas}$ for the more irradiated models ($\chi_{ISM} > 10^3 G_0$) makes the line emitting regions radially more extended, and the mid-IR line fluxes increase (Fig. A.5g).

In the outer disk, the opacity in the emitting region of the far-IR lines drops of an order of magnitude at $\chi_{ISM} = 10^4 G_0$, but then it raises again for larger $\chi_{ISM}$. This is because of a complex variation in the total water column density in the outer disk (Figs. A.5b and a), which also affects the emitting region's extension and position (Figs. A.5c and d). The far-IR lines first moves outward and then above $\chi_{ISM} = 10^3 G_0$ inward and closer to the midplane. The line optical depth drops below 1.0 for $\chi_{ISM} > 10^3 G_0$ (Fig. A.5d). Increasing the FUV radiation field, the emitting column density becomes larger because of enhanced photodesorption (Fig. A.5a). The emitting region, however, is radially less extended because of increasing photodissociation, which removes water from the outer disk and even partially depletes the ice reservoir (Fig. A.5b). The gas temperature in the emitting region increases from less than 100 K to 500 K and decouples from the dust temperature for models with $\chi_{ISM} > 100 G_0$. The interplay between the increased temperature, the larger column density and the decreasing emitting area in the outer disk limits the sub-mm line flux increase to an order of magnitude.

The theoretical *Spitzer* spectrum (Figs. A.5e and f) shows an increase of the line-to-continuum ratio for an increasing $\chi_{ISM}$ up to $10 G_0$. For fluxes larger than 1000 $G_0$, the continuum starts to become stronger, reducing the line-to-continuum ratio. This is due to the fact that $\chi_{ISM} > \chi_{star}$ ($\chi_{star} \sim 10^4 G_0$ within 2 AU). Beyond $10^3 G_0$, the mid-IR color becomes much smaller than 1.

## 4.7. PAHs photoelectric heating

PAHs affect the total water column density by less than an order of magnitude in the outermost regions of the disk (Fig. A.6a). The *Spitzer* transitions are almost unaffected (Fig. A.6b). Their emitting region moves slightly upward and outward due to the increase in heating caused by PAHs (Fig. A.6c). However, this displacement is too small to affect the line fluxes in the mid-IR. The column density in the emitting region of the far-IR lines increases by less than a factor of two with increasing $f_{PAH}$. The sub-mm line emitting region moves outward and upward (Fig. A.6d), but the gas temperature hardly changes. The line fluxes increase by a factor ~2.

Since the dust opacities stay constant, the *Spitzer* spectra all show the same continuum; the 13.5/30 µm continuum flux ratio is thus also unaffected. The line-to-continuum ratio grows with increasing PAH abundance because of the heating effect on the gas.
## 4.8. Disk gas mass

The total gas mass in the disk affects almost proportionally the total water column density in the disk (Fig. A.7b). All reservoirs are enhanced by the same amount with increasing gas mass. In the *Spitzer* emitting region, the gas density progressively increases for larger $M_{\text{gas}}$. The water column density in the *Spitzer* water line emitting region follows a nonlinear trend (Fig. A.7a). From $M_{\text{gas}} 10^{-5}$ to $10^{-4} M_{\odot}$, the emitting region moves outward and upward (Fig. A.7c). Then, for masses larger than $M_{\text{gas}} = 10^{-4} M_{\odot}$, the migration of the line emitting region is inverted. Below a certain water column density ($N_{\text{H}_2\text{O}} = 1.6-8 \times 10^{17} \text{ cm}^{-2}$, Bethell & Bergin 2009), self-shielding becomes inefficient and water is strongly photodissociated. For the *Spitzer* line fluxes, the driving mechanism is first the increase of gas column density in the emitting region and, subsequently, the growth of the emitting region itself.

The water column density in the outer disk increases with $M_{\text{gas}}$ (Fig. A.7b) and also the optical depth of the lines. This causes the transition to be emitted progressively outward and higher up in the outer disk (Fig. A.7d). Thus, the increase in the far-IR line flux with $M_{\text{gas}}$ has the same explanation as the increase in the mid-IR lines.

The theoretical *Spitzer* spectrum (Figs. A.7e and f) cleary shows an increase in the line fluxes with the gas content of the disk. The different continuum level is a consequence of the dust settling. In the formalism adopted here, the corresponding dust scale height depends on the gas density, which is directly proportional to the mass of the gas in the disk. This means that increasing $M_{\text{gas}}$ produces less settling, and hence a stronger IR continuum with a steeper slope. As a consequence, the mid-IR color becomes smaller.

# 5. Discussion

## 5.1. Discrimination between parameters affecting water line fluxes

The main goal of this work can be summarized in two questions:

1. Why do disks sometimes lack the detection of water?

2. Why are the water lines from the outer disk absent in disks that have clear detections of inner water reservoirs?

To answer these questions, we divide the parameters into three categories: those that strongly affect the line ratios (line flux variation $\ge 1$ dex, parameter variation $< 1$ dex); those that produce a relevant variation (line flux variation same order of magnitude as parameter variation); and parameters that produce a weak variation (line flux variation less than parameter variation). An overview of the effects of the parameters on continuum and line fluxes is listed in Table 7.

The *Spitzer* transitions around $12 \ \mu\text{m}$ are not those typically detected since their flux is weaker compared to other blends. We verified, however, that the same conclusions also hold for other mid-IR lines around 15, 17, and $29 \ \mu\text{m}$. The conclusions about the far-IR lines are limited to the fundamental ortho and para water lines. The line forming regions of other far-IR lines suffer very large radial and vertical displacements in these model series, hence, making the discussion even more complex. They will require a dedicated study.

The *Herschel*/HIFI lines would be largely undetectable with a standard observational setup and 1 h of exposure, however, our modeling results suggest that the lines become detectable in many disks at a $5\sigma$ sensitivity level of $5 \times 10^{-20} \text{ W/m}^2$.

<figure><img src="image_13.png" /><figcaption>**Fig. 5.** Overview of the line ratios *Herschel*/Spitzer (538.29 µm/12.407 µm) for all the parameters as a function of the explored range of values.</figcaption></figure>

TW Hya required a total 181 min on-source integration with *Herschel*/HIFI to detect the 538.29 µm line (Hogerheijde et al. 2011).

In the first category, we find parameters such as the flaring index $\beta$ and the power law dust size distribution exponent $a_{\text{pow}}$. These parameters are hard to constrain from observations of SEDs alone. The *Herschel* and *Spitzer* line fluxes are very sensitive to their exact values. Flaring affects the *Herschel*/HIFI lines strongly (2 dex) and the *Spitzer* transitions only marginally (less than a factor 3). The exponent of the dust grain size distribution produces larger changes in the *Spitzer* regime than in the *Herschel* regime. These flux variations are due to geometrical effects (in the outer disk, because of flaring, the emitting column of water and the area of the emitting region increase) and opacity/temperature changes (the inner disk in the flared models is directly exposed to the radiation, and so is warmer). In the case of the dust size distribution exponent, the far-IR is always emitted in an optically thin region, and therefore are less affected by any opacity change, while in the mid-IR regime the variation is about 4 dex and the line flux is suppressed. The external part (at larger radii) of reservoir B and all of reservoir C are optically thin in the continuum. This means that the line flux is directly affected by the variation in the dust opacity, and partially by a change in the local thermal conditions. For the transitions produced in the inner disk (as the *Spitzer* ones), the situation is driven by a change in the thermal structure since the lines are optically thick, with the exception of very peculiar models explored in our series of models.

In the second category, we find almost all other parameters: dust-to-gas mass ratio, *d/g*; maximum dust size, $a_{\text{max}}$; ISM radiation field, $\chi_{\text{ISM}}$; and gas mass of the disk, $M_{\text{gas}}$. Excluding the ISM radiation field, that is an environmental property, all the other properties are good candidates to understand the absence of a correlation between detections in the mid-IR and far-IR wavelength range. In these cases, the *Spitzer* lines are more affected by the parameter variations, and are strongly enhanced (3–4 dex) in disks with $M_{\text{dust}} < 10^{-4} M_{\odot}$ and/or $M_{\text{gas}} > 0.01 M_{\odot}$, dust grains larger than 400 µm and few small grains ($a_{\text{pow}} < 3.5$). The continuum is enhanced in the presence of larger dust masses or smaller dust grains that absorb more FUV radiation and become warmer. *Herschel* lines show the same trend, but the magnitude of the variation is smaller. The effect of these parameters
<table><caption>**Table 7.** General overview of the series of models results.</caption><thead><tr><th>Parameter*</th><th>IRS lines</th><th>mid-IR continuum</th><th>HIFI lines</th><th>IRS line profiles</th><th>HIFI line profiles</th></tr></thead><tbody><tr><td>d/g</td><td>--</td><td>+ (bolometrically)</td><td>--</td><td>narrower</td><td>unaffected</td></tr><tr><td>β</td><td>+</td><td>+ (bolometrically, far-IR stronger)</td><td>++</td><td>narrower</td><td>unaffected</td></tr><tr><td>α<sub>set</sub></td><td>--</td><td>+ (stronger with increasing λ)</td><td>--</td><td>unaffected</td><td>unaffected</td></tr><tr><td>a<sub>max</sub></td><td>++</td><td>--</td><td>+</td><td>narrower</td><td>unaffected</td></tr><tr><td>a<sub>pow</sub></td><td>--</td><td>+ (variation weaker around 15 µm)</td><td>--</td><td>wider</td><td>unaffected</td></tr><tr><td>χ<sub>ISM</sub></td><td>++<sup>1</sup></td><td>++<sup>1</sup></td><td>+</td><td>narrower</td><td>unaffected</td></tr><tr><td>f<sub>PAH</sub></td><td>x</td><td>x</td><td>+</td><td>unaffected</td><td>unaffected</td></tr><tr><td>M<sub>gas</sub></td><td>++</td><td>+ (stronger with increasing λ)</td><td>++</td><td>narrower</td><td>narrower</td></tr></tbody></table>

Notes. All the variations refer to an increase of the respective parameter value. Symbols meanings: + increase < 1 dex, – decrease < 1 dex, ++ increase ≥ 1 dex, – decrease ≥ 1 dex, x constant; (1) for $χ_{ISM} > 10^4 G_0$.

on the spectrum of the outer disk is not as important as in the inner disk.

The last category of parameters, mixing coefficient $α_{set}$ and fraction of PAHs $f_{PAH}$, do not significantly affect the *Spitzer* nor the *Herschel* water emission lines. All water transitions respond in the same manner to the mixing coefficient. We stress here that the treatment of the mixing is disconnected from the viscous heating, which is set to zero in all our models. The PAHs are more important in the outer disk, which is optically thin, and is thus be warmed up by the increased photoelectric heating from PAHs. PAHs are not treated consistently as an opacity source, otherwise the enhanced opacity (especially in the UV-optical regime) would have produced effects comparable to an increased dust-to-gas mass ratio in the disk.

Since every parameter affecting the line flux also modifies the emitting region position and extension in the disk, the line profiles are consequently affected, and outward displacements of the emitting regions produce narrower transitions (Keplerian velocity decreases with radial distances).

Our work finds that the main driver of water spectroscopy in the mid- to far-IR wavelengths range is the continuum opacity. Changes occur because of differences in the disk dust content or the opacity function ($a_{max}$, $a_{pow}$ variation). These directly affect the line flux (extinction) or indirectly affect the local thermal conditions (temperature change). The second effect is the disk geometry. A flared disk geometry exposes more water to the central stellar radiation field. Hence the total column density is depleted in the outer disk because of photodissociation. An increase in the ISM FUV radiation field also produces a suppression of water in the outermost reservoir C.

As discussed in Sect. 3.2, the mid-IR water lines are in LTE, while the outer disk water lines are marginally in non-LTE. Hence the far-IR line fluxes can be affected by uncertainties in the collision cross sections, as discussed in Kamp et al. (2013). We perform a sensitivity check on the standard model making a random variation (increase/decrease) of an order of magnitude in the water collisional constants. The results confirm the robustness of the conclusions presented in this paper: the line fluxes of the *Spitzer* transitions are unaffected, while the *Herschel* line variations are within 10%.

The disk water distribution found from the modeling is in qualitative agreement with recent modeling performed by Du & Bergin (2014). In particular, we reproduce the change of the column densities of water in the inner and outer disk due to the variation of dust-to-gas mass ratio and settling. A quantitative comparison is not possible because we use a more physical approach in our models (non-parametric settling; Dubrulle et al. 1995).

## 5.2. Comparison with observations

Observationally, it will be difficult to assess the exact properties of a certain disk. Previous observational studies of mid-IR water spectroscopy used continuum diagnostics in the same spectral region to connect to disk properties of the targets. We expect to find such correlations, and we want to use them to define the position of each model produced by a parameter space exploration (that starts from our standard disk model), in a plot of the line flux vs. continuum flux. Here, we decide to use a mid-IR color (13.5/30; similar to what is used in Acke et al. 2009) because it traces slope changes in the *Spitzer* spectrum, and we found that that every parameter produces its own fingerprint. For the far-IR, we decided to use a diagnostic connected with the outer disk properties. We choose the SCUBA flux at 850 µm because it is available for many objects and traces dust mass to first order. We search for a direct relation between mid-IR colors, submm photometry, and simulated line fluxes at 12.407 µm and 538.29 µm. We distinguished three types of parameters (not to be confused with the categories defined in Sect. 5): parameters that change the dust opacity function ($a_{max}$, $a_{pow}$), parameters that change only the total opacity ($M_{gas}$, $d/g$, $α_{set}$), and parameters that do not affect the opacity (ISM radiation field, $f_{PAH}$, $β$). Every parameter produces its own signature on the SED and line fluxes, and Fig. 6 and Table 8 summarize the main findings. The parameter $M_{gas}$ produces a limited variation in the mid-IR continuum, and a more important effect on the submm (82% in sub-mm, 10% in mid-IR, due to the settling). It produces a strong variation in the line flux both in the mid- and far-IR. Settling only modestly affects the line fluxes in both the regimes, but the mid-IR slope is strongly changed by up to 50% and the submm continuum flux changes by 87%. The dust-to-gas mass ratio does not strongly affect the mid-IR continuum (about 25%); however it changes the submm flux by almost a factor 2. Flaring significantly affects both the mid-IR continuum (90% relative variation) and the submm (more than 95% between extreme models), but the *Spitzer* lines are hardly affected. The ISM FUV radiation flux produces a stronger change in the *Spitzer* lines compared to the *Herschel*/HIFI lines, and a steepening of the IR excess beyond 10 µm. For the far-IR lines, the flux grows with increasing FUV radiation field until a certain limit and then levels off. The fraction of PAHs does not affect the SED, since the PAH opacity is not taken into account in the continuum radiative transfer.

A comparison of the modeling results with observations toward T Tauri disks shows that our simplified model predictions define a parameter space that is in good agreement with a number of objects: AA Tau, DG Tau, RW Aur, DM Tau for the mid-IR spectral region (top plot of Fig. 6), DG Tau, DM Tau,
<figure><img src="image_14.png" /><figcaption>**Fig. 6.** Line fluxes vs. continuum. Top: 12.407 µm versus 13.5/30 IRS continuum flux ratios. Bottom: 538.29 µm versus the 850 µm SCUBA flux. Color code: blue for parameters not affecting the opacity, red for parameters affecting the total continuum opacity, green for parameters affecting the opacity function. Overplotted are the sensitivity limits for *Spitzer*/IRS, *Herschel*/HIFI, and JWST/MIRI (references in the text). The magenta data points are *Spitzer*/IRS observations available for stars with spectral type similar to our model central star, and blend fluxes have been extracted from *Spitzer* spectra following the procedure described in Appendix D. Cyan data points are the only *Herschel*/HIFI observations published for this submm line. All line fluxes have been scaled to 140 pc. Arrows indicate upper limits.</figcaption></figure>

and TW Hya for the far-IR spectral region (bottom plot of Fig. 6). These targets have spectral types and properties compatible with our central star, and are not heavily affected by the presence of companions, either because they are too faint (case of RW Aur A+B; White & Ghez 2001) or because the members are widely separated (DG Tau A+B; Mundt & Fried 1983). Appendix D describes how the 12.4 µm blend flux was extracted from *Spitzer* spectra for the above targets and how continuum errors were propagated; fluxes are listed in Table D.1. Figure 6 compares the model results with the observational data: line blends around 12.4 µm versus our mid-IR color and 538.29 µm versus the SCUBA 850 µm flux. Our model series does not pretend to describe the situation of each target, which is a complex interplay of properties, but the position of the targets agree nicely with our models and support the predictive power of mid- and far-IR observations. The direct comparison with the sensitivities

of *Spitzer*/IRS and JWST/MIRI (or *Herschel*/HIFI for the far-IR plot) shows the discovery potential of these facilities.

Figure 6 shows that it is not possible to build trends that clearly connect line fluxes with continuum properties because different parameters affect the relation between these two quantities in a different manner. Hence, it is logical to expect a large amount of scatter in observations. Previous observational work (e.g., Pontoppidan et al. 2010b) show several plots with scattered observations in the continuum properties or line vs. continuum fluxes or even continuum versus continuum. Any individual parameter produces a complex displacement in each of these plots, and every disk can be a different combination of parameters.

### 5.3. The future with JWST

The low resolution of the IRS spectrograph (R = 600) does not enable an investigation of the line profiles; the blends usually contain several transitions (Banzatti 2013). The problem is particularly important when the line profiles deviate from Keplerian rotation because of other kinematical components caused by winds/outflows/jets (Baldovin-Saavedra et al. 2011, 2012). The future facility to study mid-IR water lines is the James Webb Space Telescope (Greenhouse et al. 2014). The larger mirror and higher spectral resolution of the MIRI spectrograph (D = 6.5 m, R ≈ 3000; Clampin 2011) will overcome the technical limitations of IRS. Figure 7 shows a comparison of the mid-IR spectrum of our standard model in the *Spitzer*/IRS range with the resolution of the IRS and MIRI spectrographs. The higher sensitivity and the spectral coverage (5–28.5 µm) similar to *Spitzer* IRS, will allow us to detect water lines in almost all the cases modeled in this study, within a 2.7 h exposure (10σ; Fig. 6). Figure 7 also shows that JWST/MIRI can make a decomposition of line blends into individual lines feasible, allowing, possibly for the first time, a quantitative analysis of the physical conditions in their emitting region.

## 6. Conclusions

We explore here the parameter space around a standard T Tauri disk model to identify the main drivers of mid- and far-IR water spectroscopy. We primarily analyze two transitions, a mid-IR line at 12.407 µm and a far-IR line at 538.29 µm, representative of the whole group of transitions considered. In our series of models, we find that the gas mass and dust-to-gas mass ratio simultaneously affect water lines in the *Spitzer* and *Herschel*/HIFI wavelength ranges. Disks with a low gas content or high dust content are undetectable in both regimes. In the first case, the lack of emitting gas explains the nondetections, in the second case, the continuum opacity hides the water reservoirs. These two parameters could explain the occurrence of disks without any water features.

Dust properties ($a_{max}$, $a_{pow}$) only affect the mid-IR transitions, and can explain the nondetections of *Spitzer* lines, while the *Herschel*/HIFI transitions would be detectable in case of several hours of integration. On the other hand, flaring strongly affects the far-IR regime; so the nondetection of far-IR water lines can also be due to disk geometry.

Models with strong external FUV radiation field ($\chi_{ISM} > 10^3 G_0$) show strong mid-IR lines with narrow profiles. *Herschel*/HIFI lines are less affected and increase only by a factor of a few due to the external irradiation.

Our choice of implementing the PAHs only as a heating source and not as an opacity source causes the PAHs fraction to
<table><caption>**Table 8.** Line flux and color variation for the different parameters changes.</caption><thead><tr><th>Parameter</th><th>Type<sup>a</sup></th><th>Extreme ratios<sup>b</sup> 538.29 µm</th><th>Flux variation 850 µm [%]</th><th>Extreme ratios<sup>b</sup> 12.407 µm</th><th>Flux ratio 13.5/30 variation</th></tr></thead><tbody><tr><td>a<sub>pow</sub></td><td>1</td><td>0.208</td><td>77.6</td><td>0.007</td><td>-0.453</td></tr><tr><td>a<sub>max</sub></td><td>1</td><td>3.216</td><td>85.0</td><td>24.507</td><td>0.119</td></tr><tr><td>M<sub>gas</sub></td><td>2</td><td>271.825</td><td>82.2</td><td>2600.789</td><td>-0.276</td></tr><tr><td>d/g</td><td>2</td><td>0.047</td><td>98.8</td><td>5.360</td><td>-0.253</td></tr><tr><td>α<sub>set</sub></td><td>2</td><td>0.456</td><td>87.5</td><td>0.454</td><td>-0.484</td></tr><tr><td>β</td><td>3</td><td>38.236</td><td>95.2</td><td>4.070</td><td>-0.911</td></tr><tr><td>χ<sub>ISM</sub></td><td>3</td><td>12.255</td><td>89.1</td><td>274.378</td><td>-0.532</td></tr><tr><td>f<sub>PAH</sub></td><td>3</td><td>2.362</td><td>0.18</td><td>1.127</td><td>0.003</td></tr></tbody></table>

**Notes.** (a) Type 1: parameters affecting global opacity and the opacity function; type 2: parameters affecting global opacity; type 3: parameters not affecting the opacity. (b) Defined as ratio between the maximal and minimal value if the flux increases and viceversa if the flux decreases. This gives the maximum variation of line flux because of the parameter range considered.

<figure><img src="image_15.png" /><figcaption>**Fig. 7.** Mid-IR spectra of the standard disk model, continuum subtracted, convolved, and properly rebinned with the resolution of *Spitzer*/IRS LH/SH modules (*R* = 600) and JWST/MIRI (*R* = 3000) and arbitrarily shifted by 0.05 Jy. The bottom panel shows a zoom in the range 12.35–12.39 µm, JWST spectrum is again shifted by 0.05 Jy. The range of MIRI is artificially extended up to 37 µm to compare the performances.</figcaption></figure>

only weakly affect both mid- and far-IR lines. Our prescription of the settling is based on a physical description of the effects of gas vertical mixing on dust grains, but it is not related to the heating associated with the mechanical turbulence in the disk;

the two phenomena are disconnected. In that case, the settling affects the inner disk physics very little.

We demonstrate the importance of future telescopes such as JWST in the investigation of water in disks. The higher sensitivity is capable of increasing dramatically the number of disks detected in water. The higher spectral resolution will permit us to detect unblended transitions, and to study in detail the physical conditions in the inner water reservoir.

**Acknowledgements.** The research leading to these results has received funding from the European Union Seventh Framework Programme FP7-2011 under grant agreement no 284405.

## References

Acke, B., Min, M., van den Ancker, M. E., et al. 2009, A&A, 502, L17

Adams, N., Smith, D., & Grief, D. 1978, Int. J. Mass Spectrometry Ion Phys., 26, 405

Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684

Andrews, S. M., & Williams, J. P. 2005, ApJ, 631, 1134

Andrews, S. M., & Williams, J. P. 2007, ApJ, 671, 1800

Anicich, V. G. 1993, J. Phys. Chem. Ref. Data, 22, 1469

Aresu, G., Kamp, I., Meijerink, R., et al. 2011, A&A, 526, A163

Baldovin-Saavedra, C., Audard, M., Güdel, M., et al. 2011, A&A, 528, A22

Baldovin-Saavedra, C., Audard, M., Carmona, A., et al. 2012, A&A, 543, A30

Bally, J., O'Dell, C. R., & McCaughrean, M. J. 2000, AJ, 119, 2919

Banzatti, A. 2013, Ph.D. Thesis, ETH Zurich

Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087

Beckwith, S. V. W., & Sargent, A. I. 1991, ApJ, 381, 250

Bergin, E. A., Hogerheijde, M. R., Brinch, C., et al. 2010, A&A, 521, L33

Bethell, T., & Bergin, E. 2009, Science, 326, 1675

Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11

Blum, J., & Wurm, G. 2008, ARA&A, 46, 21

Boulanger, F., Abergel, A., Bernard, J.-P., et al. 1996, A&A, 312, 256

Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859

Brogi, M., de Kok, R. J., Birkby, J. L., Schwarz, H., & Snellen, I. A. G. 2014, A&A, 565, A124

Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437

Carr, J. S., & Najita, J. R. 2008, Science, 319, 1504

Carr, J. S., & Najita, J. R. 2011, ApJ, 733, 102

Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29

Clampin, M. 2011, in SPIE Conf Ser., 8146, 814605

Crouzet, N., McCullough, P. R., Deming, D., & Madhusudhan, N. 2014, ApJ, 795, 166

de Graauw, T. 2004, BAAS, 36, 825

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133

Dent, W. R. F., Greaves, J. S., & Coulson, I. M. 2005, MNRAS, 359, 663

Di Francesco, J., Johnstone, D., Kirk, H., MacKenzie, T., & Ledwosinska, E. 2008, ApJS, 175, 277

Dominik, C., Ceccarelli, C., Hollenbach, D., & Kaufman, M. 2005, ApJ, 635, L85

Draine, B. T. 1978, ApJS, 36, 595

Draine, B. T. 2003, ARA&A, 41, 241

Draine, B. T. 2006, ApJ, 636, 1114
Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Du, F., & Bergin, E. A. 2014, ApJ, 792, 2

Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237

Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971

Dullemond, C. P., & Dominik, C. 2008, A&A, 487, 205

Encrenaz, T. 2008, ARA&A, 46, 57

Evans, II, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Fahey, D. W., Fehsenfeld, F. C., Ferguson, E. E., & Viehland, L. A. 1981, J. Chem. Phys., 75, 669

Faure, A., & Josselin, E. 2008, A&A, 492, 257

Field, D., Adams, N. G., & Smith, D. 1980, MNRAS, 192, 1

Flower, D. R., & Pineau des Forets, G. 1995, MNRAS, 275, 1049

Flower, D. R., & Pineau des Forêts, G. 2003, MNRAS, 343, 390

Fraine, J., Deming, D., Benneke, B., et al. 2014, Nature, 513, 526

Geers, V. C., van Dishoeck, E. F., Pontoppidan, K. M., et al. 2009, A&A, 495, 837

Glasse, A. C. H., Bauwens, E., Bouwman, J., et al. 2010, in SPIE Conf. Ser., 7731, 77310K

Glassgold, A. E., Meijerink, R., & Najita, J. R. 2009, ApJ, 701, 142

Green, S., Maluendes, S., & McLean, A. D. 1993, ApJS, 85, 181

Greenhouse, M. A., Kimble, R. A., Rauscher, B. J., et al. 2014, in AAS Meet. Abstr., 223, 149.29

Guilloteau, S., Di Folco, E., Dutrey, A., et al. 2013, A&A, 549, A92

Gürtler, J., Schreyer, K., Henning, T., Lemke, D., & Pfau, W. 1999, A&A, 346, 205

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Heinzeller, D., Nomura, H., Walsh, C., & Millar, T. J. 2011, ApJ, 731, 115

Herbst, E., Millar, T. J., Wlodek, S., & Bohme, D. K. 1989, A&A, 222, 205

Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338

Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2012, in Proc. Symp. From Atoms to Pebbles: Herschel's view of Star and Planet Formation, ed. J.-C. Augereau

Hollenbach, D., & Gorti, U. 2009, ApJ, 703, 1203

Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18

Isella, A., Testi, L., & Natta, A. 2006, A&A, 451, 951

Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213

Jensen, M. J., Bilodeau, R. C., Safvan, C. P., et al. 2000, ApJ, 543, 764

Kamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991

Kamp, I., Tilling, I., Woitke, P., Thi, W.-F., & Hogerheijde, M. 2010, A&A, 510, A18

Kamp, I., Woitke, P., Pinte, C., et al. 2011, A&A, 532, A85

Kamp, I., Thi, W.-F., Meeus, G., et al. 2013, A&A, 559, A24

Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714

Kim, J. K., & Huntress, Jr., W. T. 1975, J. Chem. Phys., 62, 2820

Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57

Lecar, M., Podolak, M., Sasselov, D., & Chiang, E. 2006, ApJ, 640, 1115

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Machida, R., & Abe, Y. 2010, ApJ, 716, 1252

Mallard, W. G., Westley, F., Herron, J. T., Hampson, R. F., & Frizell, D. H. 1994, NIST Chemical Kinetics Database: v. 6.0 (Gaithersburg, MD: National Institute of Standards and Technology)

Mannings, V., & Emerson, J. P. 1994, MNRAS, 267, 361

Marconi, A., Testi, L., Natta, A., & Walmsley, C. M. 1998, A&A, 330, 696

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Mauclaire, G., Derai, R., & Marx, R. 1978, Int. J. Mass Spectrometry Ion Phys., 26, 289

McClure, M. K., Manoj, P., Calvet, N., et al. 2012, ApJ, 759, L10

McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36

McKee, C. F., & Williams, J. P. 1997, ApJ, 476, 144

Meijerink, R., Poelman, D. R., Spaans, M., Tielens, A. G. G. M., & Glassgold, A. E. 2008, ApJ, 689, L57

Meijerink, R., Pontoppidan, K. M., Blake, G. A., Poelman, D. R., & Dullemond, C. P. 2009, ApJ, 704, 1471

Meijerink, R., Aresu, G., Kamp, I., et al. 2012, A&A, 547, A68

Millar, T. J., Farquhar, P. R. A., & Willacy, K. 1997, A&AS, 121, 139

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96

Mundt, R., & Fried, J. W. 1983, ApJ, 274, L83

Najita, J. R., Ádámkovics, M., & Glassgold, A. E. 2011, ApJ, 743, 147

Najita, J. R., Carr, J. S., Pontoppidan, K. M., et al. 2013, ApJ, 766, 134

Natta, A., Testi, L., Neri, R., Shepherd, D. S., & Wilner, D. J. 2004, A&A, 416, 179

Ormel, C. W., Min, M., Tielens, A. G. G. M., Dominik, C., & Paszun, D. 2011, A&A, 532, A43

Padgett, D. L., Brandner, W., Stapelfeldt, K. R., et al. 1999, AJ, 117, 1490

Piétu, V., Guilloteau, S., & Dutrey, A. 2005, A&A, 443, 945

Podio, L., Kamp, I., Codella, C., et al. 2013, ApJ, 766, L5

Pontoppidan, K. M., Salyk, C., Blake, G. A., & Käufl, H. U. 2010a, ApJ, 722, L173

Pontoppidan, K. M., Salyk, C., Blake, G. A., et al. 2010b, ApJ, 720, 887

Prasad, S. S., & Huntress, Jr., W. T. 1980, ApJS, 43, 1

Rakshit, A. 1982, Int. J. Mass Spectrometry Ion Phys., 41, 185

Riviere-Marichalar, P., Ménard, F., Thi, W. F., et al. 2012, A&A, 538, L3

Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, J. Quant. Spectr. Rad. Transf., 130, 4

Salyk, C., Pontoppidan, K. M., Blake, G. A., et al. 2008, ApJ, 676, L49

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Sicilia-Aguilar, A., Henning, T., Dullemond, C. P., et al. 2011, ApJ, 742, 39

Smith, D., Panel, P., & Mayhew, C. A. 1992, Int. J. Mass Spectrometry Ion Processes, 117, 457

Smith, N., Bally, J., Licht, D., & Walawender, J. 2005, AJ, 129, 382

Stapelfeldt, K. R., Burrows, C. J., Krist, J. E., et al. 1998, ApJ, 508, 736

Tennyson, J., Zobov, N. F., Williamson, R., Polyansky, O. L., & Bernath, P. F. 2001, J. Phys. Chem. Ref. Data, 30, 735

Terada, H., Tokunaga, A. T., Kobayashi, N., et al. 2007, ApJ, 667, 303

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI (Tucson: University of Arizona Press), 339

Tielens, A. G. G. M. 2008, ARA&A, 46, 289

Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722

van der Tak, F. F. S., Chavarría, L., Herpin, F., et al. 2013, A&A, 554, A83

van der Werf, P. P., Isaak, K. G., Meijerink, R., et al. 2010, A&A, 518, L42

van Dishoeck, E. F., Jonkheid, B., & van Hemert, M. C. 2006, Faraday Discussions, 133, 231

van Dishoeck, E. F., Bergin, E. A., Lis, D. C., & Lunine, J. I. 2014, Protostars and Planets VI, 835

Weidenschilling, S. J. 1997, Icarus, 127, 290

White, R. J., & Ghez, A. M. 2001, ApJ, 556, 265

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Wilner, D. J., D'Alessio, P., Calvet, N., Claussen, M. J., & Hartmann, L. 2005, ApJ, 626, L109

Woitke, P. 2012, in Proc. Symp. From Atoms to Pebbles: Herschel's view of Star and Planet Formation, ed. J.-C. Augereau, 39

Woitke, P., Kamp, I., & Thi, W.-F. 2009a, A&A, 501, 383

Woitke, P., Thi, W.-F., Kamp, I., & Hogerheijde, M. R. 2009b, A&A, 501, L5

Woitke, P., Riaz, B., Duchêne, G., et al. 2011, A&A, 534, A44

Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197

Wyatt, M. C. 2008, ARA&A, 46, 339

Zhang, K., Pontoppidan, K. M., Salyk, C., & Blake, G. A. 2013, ApJ, 766, 82

Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481

Zsom, A., Sándor, Z., & Dullemond, C. P. 2011, A&A, 527, A10

Zuckerman, B., Forveille, T., & Kastner, J. H. 1995, Nature, 373, 494
Appendix A: Additional figures

<figure><img src="image_16.png" /><figcaption>Fig. A.1. Mixing coefficient models. For further details see Fig. 4.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. A.2. Flaring index exponent models. For further details see Fig. 4.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Fig. A.3. Maximum dust size models. For further details see Fig. 4.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. A.4. Dust size distribution power law index models. For further details, see Fig. 4.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. A.5. ISM radiation field models. For further details see Fig. 4.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Fig. A.6.** Fraction of PAHs models. **a)** Emitting region water column density for dust-to-gas mass ratio models; this region comprises 15–85% of the radially and vertically integrated flux. Green vertical dotted lines distinguish the reservoirs. **b)** Theoretical *Spitzer* SH/LH modules spectra (R = 600) for a dust-to-gas model series, continuum subtracted, and arbitrarily shifted. **c)** Average extension of the line emitting region, which comprises 50% of the line flux radially and vertically for 12.407 µm. **d)** Average extension of the line emitting region, which comprises 50% of the line flux radially and vertically for 538.29 µm.</figcaption></figure>
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /><figcaption>Fig. A.7. Gas mass content models. For further details see Fig. 4.</figcaption></figure>
## Appendix B: Computation of properties in the line emitting region

The different water lines are divided in groups based on their spatial origin. The latter is defined as that part of the disk from which a cumulative radially and vertically integrated flux of 15% to 85% of the total flux is produced. For each of these regions, averaged physical conditions $\langle X \rangle$ are derived, such as the density of water, density of the collisional partners (e-, H, H$_2$), $T_{gas}$ and $T_{dust}$ by numerically integrating over the emitting region (Eq. (B.1)), i.e.,

$$ \langle X \rangle = \frac{\int_{r_{\text{in}}}^{r_{\text{out}}} \int_{z_{\text{low}}}^{z_{\text{high}}} X n_{\text{gas}} 2\pi r dr dz}{\int_{r_{\text{in}}}^{r_{\text{out}}} \int_{z_{\text{low}}}^{z_{\text{high}}} n_{\text{gas}} 2\pi r dr dz} \quad (B.1) $$

In order to determine the continuum extinction, a vertical integration is performed from the disk surface to the midplane. $\tau_v$ is the continuum optical depth at the frequency of the line. It is related to the dust density $\rho_{dust}$ and dust mass extinction coefficient $\kappa_v$. The integration is performed numerically, as follows:

$$ \tau_v(z) = \int_{z_{\text{max}}(r)}^{z(r, \tau_v=1)} \rho(r, z) \kappa_v dz. \quad (B.2) $$

Through interpolation, we find the value at which $\tau_v = 1$. The optical depth is evaluated at the radial midpoint of the emitting region.

## Appendix C: Model details

The following appendix shows a number of additional plots. These can be useful for the interpretation of the main results and conclusions reported in the paper.

Here we show the Herschel/HIFI fundamental o-H$_2$O and p-H$_2$O lines fluxes behavior with the different parameters.

Here, we also present plots collecting the main properties of the standard model: $T_{gas}$, $T_{dust}$, the local dust-to-gas mass ratio, local average dust size, and the main heating and cooling processes.

Line ratio plots are used to compare the relative behavior of the mid- and far-IR lines, considering two representative Spitzer water transitions and the Herschel fundamental ortho and para water lines. To show the trend clearly, we add a parabolic regression to the data (magenta lines).
<figure><img src="image_37.png" /></figure>

**Fig. C.1.** Line fluxes behavior for Herschel/HIFI lines 538.29 µm and 269.27 µm for different model series. First row left: dust-to-gas mass ratio, d/g. First row right: mixing parameter, α$_{set}$. Second row left: disk flaring, β. Second row right: maximum grain size a$_{max}$. Third row left: power law index of grain size distribution, a$_{pow}$. Third row right: strength of UV radiation field, χ$_{ISM}$. Bottom left: PAH abundance, f$_{PAH}$. Bottom right: disk gas mass, M$_{gas}$.
<figure><img src="image_38.png" /><figcaption>**Fig. C.2.** Standard model, from top to bottom and from left to right: *T*<sub>gas</sub> with overplotted gas temperature contours (white/black) and *A*<sub>V</sub> = 1, 10 (red), *T*<sub>dust</sub> with overplotted temperature contours (black/white) and *A*<sub>V</sub> = 1, 10 (red), dust-to-gas mass ratio, average dust grain size, main cooling processes (CII fine structure lines, dust grain-gas thermal exchange, Ly α line, neutral O line cooling, CO rotational and ro-vibrational cooling, water rotational and ro-vibrational cooling, formaldehyde rotational cooling, ammonia rotational cooling, OIII line cooling), main heating processes (collisional excitation of vibrationally excited H<sub>2</sub>, chemical heating due to H<sub>2</sub> formation on dust grains, dust grain-gas thermal exchange, heating from the superthermal photo-electrons produced by C, heating by cosmic rays, heating by the photo-electrons produced by polycyclic aromatic hydrocarbons, IR heating due to the water transitions, heating by mutual electron collisions, Coulomb heating, produced by X-rays photoionization, reaction free energy heating from exothermic reactions).</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>**Fig. C.3.** Ratios between Herschel/HIFI & Spitzer/IRS lines 538.29 µm/12.407 µm and 269.27 µm/12.453 µm for different model series. First row: dust-to-gas mass ratio, *d/g*. Second row: mixing parameter, *α*<sub>set</sub>. Third row: disk flaring, *β*. Fourth row: maximum grain size, *a*<sub>max</sub>.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>**Fig. C.4.** Ratios between Herschel/HIFI & Spitzer/IRS lines 538.29 µm/12.407 µm and 269.27 µm/12.453 µm for different model series. Top: power law index of grain size distribution, *aPow*. Second line: strength of UV radiation field, $\chi_{ISM}$. Third line: PAH abundance, $f_{PAH}$. Fourth line: disk gas mass, $M_{gas}$.</figcaption></figure>
## Appendix D: Mid- and far-IR observations

Far-IR line fluxes have been obtained from the literature. Mid-IR blends have been extracted from IRS spectra fitting a Gaussian function plus a linear continuum, computing the error in the flux as noise of the continuum,

$$F(\lambda) = A \times e^{(\lambda-\lambda_0)^2/(2\sigma)^2} + B \times \lambda + C. \quad (D.1)$$

This procedure is formally identical to that used in Pontoppidan et al. (2010b). Table D.1 lists all the observational data shown in Fig. 6 with respective references.

The mid-IR color 13.5/30 continuum flux is defined as pure ratio between the Spitzer/IRS fluxes at 13.5 µm and 30 µm, similarly to what is done in Acke et al. (2009) for the 30/13.5 continuum ratio. The observations overplotted include the errors propagated accordingly to

$$\sigma_{\text{ratio}} = \frac{F_{13.5}}{F_{30}} \cdot \sqrt{\left(\frac{\sigma_{13.5}}{F_{13.5}}\right)^2 + \left(\frac{\sigma_{30}}{F_{30}}\right)^2}. \quad (D.2)$$

<table><caption>**Table D.1.** Line and continuum observations of the targets discussed in the main paper.</caption><thead><tr><th>Object</th><th><i>d</i> [pc]</th><th>12.407 µm blend flux<sup>**</sup><br>[10<sup>-17</sup> W/m<sup>2</sup>]</th><th>13.5/30 µm continuum<br>flux ratio<sup>1</sup></th><th>538.29 µm line flux<br>[10<sup>-20</sup> W/m<sup>2</sup>]</th><th>850 µm continuum<sup>5</sup><br>[Jy]</th></tr></thead><tbody><tr><td>AA Tau</td><td>140</td><td>&lt;1.33</td><td>0.663 ± 0.027</td><td>n.a.</td><td>n.a.</td></tr><tr><td>DM Tau</td><td>140</td><td>n.a.</td><td>n.a.</td><td>&lt;1.15<sup>4</sup></td><td>0.237 ± 0.012</td></tr><tr><td>DG Tau</td><td>140</td><td>&lt;10.63</td><td>0.388 ± 0.008</td><td>67.00 ± 7.00<sup>2</sup></td><td>1.180 ± 0.118</td></tr><tr><td>RW Aur</td><td>144</td><td>&lt;13.51<sup>*</sup></td><td>0.702 ± 0.010<sup>*</sup></td><td>n.a.</td><td>n.a.</td></tr><tr><td>TW Hya</td><td>56</td><td>n.a.</td><td>n.a.</td><td>2.72 ± 0.12<sup>3, *</sup></td><td>0.182 ± 0.018<sup>*</sup></td></tr></tbody></table>

**Notes.** (1) Ratios have been computed from c2d IRS spectra (Evans et al. 2009); (2) from Podio et al. (2013); (3) scaled to 140 pc from Hogerheijde et al. (2011); (4) from Bergin et al. (2010); (5) scaled to 140 pc from Di Francesco et al. (2008); (*) flux scaled to 140 pc; (**) 12.407 µm fluxes extraction procedure is described in the text of this Appendix.