WATER VAPOR DISTRIBUTION IN PROTOPLANETARY DISKS

The layout of our code¹ is similar to the ProDiMo code (Woitke et al. 2009a). Given a distribution of gas and dust, we first solve the dust temperature distribution with a Monte Carlo method based on the strategy of Bjorkman & Wood (2001; see also Baes et al. 2005; Bruderer et al. 2012), in which the dust temperature of a spatial cell is updated each time a photon packet cross this cell; this strategy is also adopted in the RADMC code (Dullemond & Dominik 2004a). Although the geometry in our code is symmetric with respect to rotation around the central axis and reflection about the midplane, the photon propagation is done in full three dimensions. We have not yet implemented the diffusion approximation (Min et al. 2009) for the highly shadowed region where the photon statistics is low, and we mitigate this by allowing a large number (∼10⁷) of photon packets in the Monte Carlo. The whole spectrum (from UV to submillimeter) of the central star is used as input. Observation and modeling of H₂ fluorescence have shown that Lyα emission can dominate the UV spectrum (Bergin et al. 2003; Herczeg et al. 2004; Schindhelm et al. 2012). In this case, the resonant scattering by atomic hydrogen in the photodissociated regions is also important, and we treated this similar to Bethell & Bergin (2011b). To increase the signal-to-noise ratio of line features like Lyα, a smaller energy is used for photon packets when its frequency falls into the line profile. We also include the absorption of UV photons by water. Since the abundance of atomic hydrogen and water is affected by chemistry, the whole process has to be iterated, which is slow but affordable. A byproduct of the radiative transfer is the distribution of radiation field over the whole disk, which will be used as input for chemistry and gas thermal balance.

After establishing the dust temperature distribution, we evolve the disk chemistry for 1 Myr. Since the heating and cooling processes are coupled with chemistry, the gas temperature is evolved in tandem with chemistry based on the heating and cooling rates. Namely, we solve the following set of ordinary differential equations (ODEs):

$$\begin{eqnarray} \frac{d}{dt} X_i &= P_i(X;T) - D_i(X;T), \quad i=1,\ldots, N, \nonumber \\ \frac{d}{dt} T &= (\Gamma - \Lambda) / C_{\rm v}, \end{eqnarray} \tag{ 1 }$$

where X_i is the abundance of species i, P_i and D_i are the production and destruction rates of this species, which are functions of the chemical abundances and temperature (and other physical parameters), and N is the total number of species. C_v = 3k_B/2 is the volume-specific heat capacity of an ideal gas, where k_B is the Boltzmann constant. The exact value of C_v is not important because we are only concerned with the equilibrium temperature rather than the rate of temperature change. The heating and cooling rates are contained in Γ and Λ. We do not need a separate set of equations to account for the elemental conservation since elements are automatically conserved within numerical tolerance. For solving the above set of ODEs, we use the DLSODES solver of the ODEPACK² package (Hindmarsh 1983), which makes use of the sparse structure of the chemical network.

The initial chemical composition is listed in Table 1. The gas temperature is set to the dust temperature at t = 0 and usually reaches steady state within a short period. We note that the chemistry cannot always reach a steady state within 1 Myr and may still evolve at a timescale ≳ 10⁸ yr.

We could also solve the chemical equilibrium (or rather quasi-equilibrium) and thermal equilibrium independently. However, in such an approach iteration for each single grid point will be needed to achieve a joint convergence, which may pose some numerical issues and takes more CPU time, while in our approach thermal equilibrium is guaranteed as far as the heating/cooling timescale is shorter that the timescale of interest (∼1 Myr).

One note about the global iteration in our code: for the radiative transfer in the first iteration, only dust is assumed to be present. This gives a distribution of dust temperature and radiation intensity over the whole disk. Based on this, the chemistry and gas temperature is solved on the grid points in a downward (i.e., from surface to midplane) then outward (from close to the central star to the disk outer edge) order. This order has the advantage that the self-shielding of H₂ and CO can be updated each time a grid point has been calculated during one iteration. The radiative transfer is redone before each chemical and thermal calculation of the whole disk to take into account the effects (H scattering, water absorption) due to updated chemical composition. The changes in the radial water abundance profile at different vertical height as the iteration proceeds will be described later in Section 3.6. Since the code has a Monte Carlo component (for the radiative transfer) built in, perfect convergence is not expected.

We use the full UMIST RATE06 network (Woodall et al. 2007) for our gas phase chemistry. Details for the implementation of this network can be seen in that paper. In addition, we include dissociation of H₂O and OH by Lyα photons, adsorption of major gas phase species onto the dust grain, and desorption of species on the dust grain surface either thermally or induced by cosmic rays and UV photons. Two-body reactions on the dust grain surface are also included, leading to the formation of H₂O, CH₃OH, CH₄, etc. The surface network is taken from Hasegawa et al. (1992). Recombination of ions with charged dust grains is included. In total, the chemical network has 467 species and 4801 reactions. We describe some reaction types of special importance in the following.

2.2.1. Adsorption

The adsorption rate of species X is

$$\begin{equation} R_{\rm ad}({\rm X}) = s \sigma v_{\rm T} n_{\rm dust}, \end{equation} \tag{ 2 }$$

where s is the sticking coefficient, σ is the cross section of dust particles, v_T is the thermal speed of species X, and n_dust is the density of dust particles. We have

$$\begin{equation} \sigma = \pi a^2,\quad v_{\rm T} = \sqrt{\frac{8k_{\rm B} T_{\rm gas}}{\pi m_{\rm X}}}, \end{equation} \tag{ 3 }$$

where a is the average radius of dust particles, k_B is the Boltzmann constant, and m_X is the mass of a particle of X.

We calculate the sticking coefficient using a formula from Chaabouni et al. (2012):

$$\begin{equation} s = \frac{1 + T_{\rm gas}/(15\tilde{m}_{\rm X})}{[1+T_{\rm gas}/(38.5\tilde{m}_{\rm X})]^{2.5}}, \end{equation} \tag{ 4 }$$

where $\tilde{m}_{\rm X}$ is the mass number of X. The numbers in the above formula are interpolated from the parameters for nonporous amorphous solid water ice and silicate dust in Chaabouni et al. (2012). The effect of using this formula rather than the commonly used constant value of one is most important for atomic H at high temperatures. If a constant value is used, the formation of H₂ may heat the gas to unrealistically high temperatures in the photodissociated region.

2.2.2. Thermal Desorption

The thermal desorption rate is

$$\begin{equation} k_{\rm evap,th} = \nu e^{-E_{\rm des}/T_{\rm dust}}, \end{equation} \tag{ 5 }$$

where ν is the characteristic vibrational frequency of species X on the dust grain surface (Hasegawa et al. 1992),

$$\begin{equation*} \nu = \sqrt{\frac{2 n_{\rm S} E_{\rm des}}{\pi ^2 m_{\rm X}}}, \end{equation*}$$

n_S being the number density of surface sites, usually taken to be 10¹⁵ cm⁻², and E_des is the desorption energy of species X, for which we adopt the values from Garrod et al. (2008). Typically ν is of the order of 10¹² Hz.

2.2.3. Photodesorption

The yield of a species on the dust grain surface per incident UV photon can be empirically written as (Öberg et al. 2009a, 2009b)

$$\begin{equation} Y = (a + b\times T_{\rm dust}) (1-e^{-x/l}), \end{equation} \tag{ 6 }$$

where x is the thickness of the ice, and a, b, and l may be approximated as constants, which are determined experimentally. The UV flux is calculated in the radiative transfer part of the code.

For H₂O, CO, and CO₂, we use the measured value of the (a, b, l) parameters from Öberg et al. (2009a, 2009b):

$$\begin{eqnarray*} \mathrm{H_2O}{\rm :}\ & 1.3\times 10^{-3}, \; 3.2\times 10^{-5}, \; 2,\\ \mathrm{CO}{\rm :} \ & 2.7\times 10^{-3}, \; 0, \; 1,\\ \mathrm{CO_2}{\rm :}\ & 2\times 10^{-3}, \; 0, \; 3. \end{eqnarray*}$$

For all the other species, we assume a = 10⁻⁴, b = 0, and l = 1.

2.2.4. Cosmic-Ray Desorption

Cosmic-ray desorption is important for the deep and cold region of the disk. We adopt the treatment of Hasegawa & Herbst (1993), in which a dust grain is assumed to be episodically heated by cosmic rays to a high temperature (70 K) followed by evaporation of a large fraction of its ice mantle. Namely,

$$\begin{equation} k_{\rm evap,CR} = k_{\rm evap,th}(T{=}70\;{\rm K})\; f(70\;{\rm K}), \end{equation} \tag{ 7 }$$

where f(70 K) is the fraction of time for the dust to spend at temperature ∼70 K, which is estimated to be 3.16 × 10⁻¹⁹ for a dust grain size of 0.1 μm and a total cosmic-ray ionization rate of ∼10⁻¹⁷ s⁻¹ (Leger et al. 1985). The cosmic-ray intensity is attenuated with an e-fold column density of 96 g cm⁻² (Umebayashi & Nakano 1981), and the induced evaporation rate is scaled down accordingly. A different grain size distribution and cosmic-ray spectrum would produce different values for these parameters, though many details are subject to large uncertainties (Cleeves et al. 2013).

2.2.5. H₂ Formation

Since H participates in many surface reactions other than the formation of H₂, we cannot simply assume that all the H atoms adsorbed onto the dust grain are converted into H₂ molecule. So we treat the formation of H₂ on dust grain surface as a normal two-body reaction between two H atoms. Hence the formation rate (i.e., the number of H₂ molecules formed per unit volume per unit time) is

$$\begin{equation} R(\mathrm{H_2}) = k_{\rm HH} n^2({\rm H}_{\rm gr}), \end{equation} \tag{ 8 }$$

where n(H_gr) is the number density of H atoms on the grain surface, which recombine with a rate coefficient

$$\begin{equation} k_{\rm HH} = \frac{\nu _{\rm H}}{n_{\rm dust}N_{\rm S}} e^{-E_{\rm diff}/T_{\rm dust}}, \end{equation} \tag{ 9 }$$

in which E_diff is the energy barrier for migrating over the dust grain surface, usually taken to be half of the desorption energy, N_S is the number of sites per dust grain, and n_dust is the number density of dust particles. At low temperatures, quantum tunneling becomes important and the exponential part will be replaced by (Hasegawa et al. 1992)

$$\begin{equation*} e^{-{2a_{\rm diff}}\sqrt{2m_{\rm X}E_{\rm diff}}/{\hbar }}, \end{equation*}$$

where a_diff of the order of 1 Å is the barrier width for surface migration. We assume H atoms on dust grain surface are chemisorbed, and set the desorption energy to 10⁴ K according to Cazaux & Tielens (2004). Physisorption alone is not enough to account for the abundance of H₂ in the hot regions of the interstellar medium (ISM; Cazaux et al. 2005). Note that if we assume all the adsorbed H atoms are converted into H₂, then R(H₂) simply becomes half of the adsorption rate of H as calculated from Equation (2).

2.2.6. Photodissociation of H₂O and OH

2.2.7. Photodissociation of CO and H₂

2.2.8. Photodissociation of Other Species

For other species in the UMIST network, their photodissociation rates are calculated based on the formula given in Woodall et al. (2007):

$$\begin{equation} G_0 \alpha \exp (-\gamma A_{\rm V}), \end{equation} \tag{ 10 }$$

where G₀ is the unattenuated UV continuum intensity at each location of the disk relative to the standard ISM value. Namely, G₀ is calculated from the stellar spectrum assuming only inverse-square-law dilution. The attenuation due to dust and possibly water is included in the A_V parameter, which is calculated by comparing the local actual UV field (obtained from Monte Carlo radiative transfer) with the unattenuated one,

$$\begin{equation} A_{\rm V} = \max (0,\;-1.086 \ln [F_{\rm UV,att}/F_{\rm UV,unatt}]). \end{equation} \tag{ 11 }$$

Such a treatment is similar to Fogel et al. (2011).

2.3.1. Photoelectric Heating

For the heating rate due to small grains and polycyclic aromatic hydrocarbons (PAHs) we use the formula from Bakes & Tielens (1994),

$$\begin{equation} \Gamma _{\rm pe} = 10^{-24} \epsilon G_0 n_{\rm H}, \end{equation} \tag{ 12 }$$

with unit erg s⁻¹ cm⁻³. is given by

$$\begin{eqnarray} \epsilon &= \frac{4.87\times 10^{-2}}{1 + 4\times 10^{-3}(G_0T^{1/2}/n_{\rm e})^{0.73}} \nonumber \\ &\quad + \frac{3.65\times 10^{-2} (T/10^4)^{0.7}}{1 + 2\times 10^{-4} (G_0 T^{1/2}/n_{\rm e})}. \end{eqnarray} \tag{ 13 }$$

The PAH abundance used by Bakes & Tielens (1994) is 1.6 × 10⁻⁷ relative to H. We take into account the effect of dust settling and growth on this heating rate by scaling down the above rate with a factor equal to the dust-to-gas mass ratio relative to the ISM value (0.01), though the actual amount of PAH in disks is uncertain.

2.3.2. Chemical Heating and Cooling

Chemical reactions can release or absorb energy. The heating/cooling due to chemical reactions can be important for keeping the model self-consistent and has been considered in some of the previous works (see, e.g., Glassgold & Langer 1973; Dalgarno & Oppenheimer 1974; Hollenbach & McKee 1979; Glassgold et al. 2012). We include the contribution to energy balance from reactions involving the major abundant species. The exothermicity or endothermicity of these reactions are calculated based on the enthalpy of the formation of the reactants and products using the following formula:

$$\begin{equation} \Delta H = \sum \nu _i\Delta _{\rm f}H^{\rm o}(i), \end{equation} \tag{ 14 }$$

where the sum is over the reactants and products in a reaction, ν_i is the stoichiometric coefficient (negative for reactants and positive for products), and Δ_fH^o(i) is the enthalpy of formation of a species. ΔH > 0 means the reaction is endothermic. The contribution of a reaction to the heating/cooling rate is kΔH, where k is the reaction rate. The enthalpy of formation of chemical species are slowly changing functions of temperature and pressure. For our purpose, it suffices to use the values measured at standard condition (i.e., p = 1 bar, T = 298 K). The thermochemical data are taken from the NIST Web book,³Binnewies & Milke (2002), Vandooren et al. (1991), and Nagy et al. (2010). In total, 591 reactions are included to contribute to the gas heating/cooling.

2.3.3. Heating by H₂ Formation

Similar to Sternberg & Dalgarno (1989) and Röllig et al. (2006), we assume one third of the energy released (=4.5 eV) in the formation of a H₂ molecule from combination of two H atoms are converted into heat of the gas. The corresponding heating rate is

$$\begin{equation} \Gamma _{{\rm H}_2\;{\rm form}} = 2.4\times 10^{-12} {\rm erg} \,R({\rm H}_2), \end{equation} \tag{ 15 }$$

where R(H₂) (see Equation (8)) is the formation rate of H₂ in unit of cm⁻³ s⁻¹.

2.3.4. Heating by Viscous Dissipation

We use the usual α prescription. The heating rate is

$$\begin{equation} \Gamma _{\rm vis} = \frac{9}{4} \alpha \rho c_{\rm S}^2 \omega _{\rm K}, \end{equation} \tag{ 16 }$$

where ρ is the mass density of the gas, c_S is the sound speed, and ω_K is the Keplerian angular velocity.

Usually α is assumed to be a constant of the order of 0.01–1. As noted by Woitke et al. (2009a), the heating rate calculated from Equation (16) can become unphysical and gives very high temperature (>10⁴ K) when the density is very low. Hence the calculated high temperature in the top layers of the disk may not be trusted, though this does not affect our goal of study, which is focused on the deeper shielded region. Thus we use the analytical formula of Bai & Stone (2011) fitted from non-ideal magnetohydrodynamical simulations:

$$\begin{eqnarray} \alpha &=& \frac{1/2}{[(2500/Am^{2.4} + (8/Am^{0.3} + 1)^2]^{1/2}}, \nonumber \\ Am &=& n_{\rm ion} \beta _{\rm ion} / \Omega, \end{eqnarray} \tag{ 17 }$$

where n_ion is the ion density, β_ion is the ion–neutral collision rate, and Ω is the local Kepler angular velocity. At the surface of the disk where the density is low, the ambipolar diffusion parameter Am will be small and so will α (∼10⁻⁴), which will partially alleviate the problem of temperature getting unphysically high.

2.3.5. Heating by Cosmic Ray and X-Ray

The cosmic-ray heating rate is (Bruderer et al. 2009a)

$$\begin{equation} \Gamma _{{\rm CR}} = 1.5\times 10^{-11} \times \zeta _{\rm CR} n_{\rm gas}, \end{equation} \tag{ 18 }$$

where ζ_CR is the cosmic-ray flux and n_gas is the gas density.

For the X-ray heating, we calculate the X-ray photoelectric cross sections for the gas and dust using the interpolation table in Bethell & Bergin (2011a) assuming a representative X-ray photon energy of 1 keV (corresponding to a 10⁷ K blackbody), and assume that each ion pair release 18 eV into the gas (Glassgold et al. 2012). The total X-ray flux from the central star is taken to be 10⁻³ L_☉. The X-ray intensity at each location is attenuated by the column toward the central star, similar to Glassgold et al. (2004; see also Glassgold et al. 1997). The contribution of X-ray to the ionization rates are treated similar to Bruderer et al. (2009b), namely, we enhance the cosmic-ray ionization rates by the calculated X-ray ionization rates. As remarked by Bethell & Bergin (2011a), the scattered diffusion of X-ray with energy ∼1–2 keV is not important, though a scattered hard X-ray field can be important even deep in the midplane, especially when the cosmic-ray intensity has been suppressed by the stellar wind of the central star (Cleeves et al. 2013). Hence our simple treatment tends to underestimate the X-ray ionization rate. Aresu et al. (2011) find that a high X-ray luminosity of the central star can reduce the amount of hot water vapor by a factor of 20 relative to the case with zero X-ray luminosity. In our test runs, we did not see such an effect, and we interpret this as (see also Meijerink et al. 2012) due to the fact that our parameterized model does not adjust the vertical structure accordingly when the heating rate is increased, which would otherwise decrease the opacity to the dissociating UV photons.

2.3.6. Energy Exchange by Gas–Dust Collision

The energy exchange between gas and dust particles due to collisions can heat or cool the gas, depending on whether the gas temperature is lower or higher than the dust. The energy exchange per unit volume per unit time is (Hollenbach & McKee 1979)

$$\begin{equation} \Lambda _{\rm gd} = 2k_{\rm B} (T_{\rm gas} - T_{\rm dust}) v_{\rm T} \sigma _{\rm d} n_{\rm dust} n_{\rm gas} f_{\rm a}, \end{equation} \tag{ 19 }$$

where the factor two is due to the fact that more energetic particles collide with the dust grain more frequently, v_T is the thermal speed of gas particles (Equation (3)), σ_d is the mean cross section of dust particles, and f_a is the accommodation coefficient, for which we take the expression from Hollenbach & McKee (1989),

$$\begin{equation} f_{\rm a} = 1 - 0.8 e^{-75/T_{\rm gas}}. \end{equation} \tag{ 20 }$$

If more than one dust species exists, their contributions are added together. We take into account the possibility that one type of dust particle heats the gas while another cools the gas.

2.3.7. Other Heating and Cooling Mechanisms

Heating by photodissociation of H₂. We use the formula from Tielens (2005):

$$\begin{equation} \Lambda _{{\rm ph},\mathrm{H_2}} = 1.36\times 10^{-23} G_0 n(\mathrm{H_2}), \end{equation} \tag{ 21 }$$

where the unit is erg s⁻¹ cm⁻³. G₀ is the local UV intensity with self-shielding and dust extinction taken into account.

Heating by photodissociation of H₂O and OH. Their contributions to heating are calculated with

$$\begin{eqnarray} \Gamma _{{\rm ph},\mathrm{H_2O}} &= F_{{\rm Ly}\alpha } \sigma _{\mathrm{H_2O}} n(\mathrm{H_2O}) E_{\mathrm{H_2O}}, \nonumber \\ \Gamma _{{\rm ph},\mathrm{OH}} &= F_{{\rm Ly}\alpha } \sigma _{\mathrm{OH}} n(\mathrm{OH}) E_{\mathrm{OH}}, \end{eqnarray} \tag{ 22 }$$

where F_Lyα is the local Lyα number flux, σ is the photodissociation cross section, and we take $E_\mathrm{H_2O}=8.1\times 10^{-12}$ erg and E_OH = 9.2 × 10⁻¹² erg, obtained by subtracting from the Lyα photon energy the enthalpy change of the two dissociation reactions. E_H2O may overestimate the actual value by a factor of two to four, since a portion of the absorbed energy can be used to excite the internal modes of OH (Mordaunt et al. 1994).

Heating by ionization of atomic carbon. We use the formula from Tielens (2005):

$$\begin{equation} \Gamma _{\rm ion,C} = 2.2\times 10^{-22} G_0 n(\mathrm{C}), \end{equation} \tag{ 23 }$$

with unit erg s⁻¹ cm⁻³.

Cooling by electrons recombine with small dust grains. We use the analytical formula from (Bakes & Tielens 1994):

$$\begin{eqnarray} \Lambda _{\rm recom} &=& 3.49\times 10^{-30}\ n_{\rm e} n_{\rm H}\ T^{0.944} \nonumber\\ &&\quad \times \left(\frac{G_0 T^{1/2}}{n_{\rm e}}\right)^{\beta }, \end{eqnarray} \tag{ 24 }$$

with unit erg s⁻¹ cm⁻³. β = 0.735/T^0.068. This cooling rate is reduced when dust is depleted.

Cooling by the rotational transitions of H₂, and the rotational and vibrational transitions of CO and H₂O. We calculate these cooling rates using the interpolation tables from Neufeld & Kaufman (1993) and Neufeld et al. (1995) rather than direct radiative transfer for the sake of speed. These tables were made based on the escape probability approximation, in which a velocity gradient is needed, for which we use the radial gradient of the orbital speed

$$\begin{equation*} {dv}/{dr} = \frac{1}{2}\left(\frac{G M}{r^3}\right)^{1/2}, \end{equation*}$$

though note that since the orbital velocity is perpendicular to the radial direction, the photons cannot easily escape in the radial direction, but rather at an angle with it.

Heating and cooling by the vibrational transitions of H₂, and cooling by C + and O emission. We use the analytical formulae from Röllig et al. (2006). The escape probability for the C + and O lines are also calculated based on the radial gradient of the orbital speed and the local velocity dispersion.

Cooling by Lyα emission, free–bound, and free–free emissions. They are usually not important for our purpose. We include them for completeness, using the formulae in Tielens (2005) and Draine (2011).

We assume that the disk is static in the vertical direction and is in Keplerian rotation in the azimuthal direction (only needed for the line radiative transfer). The axisymmetric disk density structure we use takes the following parameterized form in cylindrical coordinates (r, z) (Lynden-Bell & Pringle 1974; Hartmann et al. 1998; Andrews et al. 2009; Cleeves et al. 2013):

$$\begin{equation} \rho (r, z) = \frac{\Sigma }{\sqrt{2\pi }h} \exp \left[-\frac{1}{2}\left(\frac{z}{h}\right)^2\right], \end{equation} \tag{ 25 }$$

where

$$\begin{eqnarray} \Sigma (r) &= \Sigma _{\rm c} \left(\frac{r}{r_{\rm c}}\right)^{-\gamma } \exp \left[-\left(\frac{r}{r_{\rm c}}\right)^{2-\gamma }\right], \nonumber \\ h &= h_{\rm c} \left(\frac{r}{r_{\rm c}}\right)^\psi. \end{eqnarray} \tag{ 26 }$$

The disk mass (gas or dust) is

$$\begin{eqnarray} M_{\rm disk} &= \int _{r_{\rm in}}^{r_{\rm out}} \Sigma \; 2\pi r dr \nonumber \\ &=\frac{2}{2-\gamma } \pi r_{\rm c}^2 \Sigma _{\rm c} \left[e^{-\left(\frac{r_{\rm in}\vphantom{_g}}{r_{\rm c}}\right)^{2-\gamma }} - e^{-\left(\frac{r_{\rm out}\vphantom{_g}}{r_{\rm c}}\right)^{2-\gamma }}\right]. \end{eqnarray} \tag{ 27 }$$

A list of the parameters involved and their meanings are in Table 2. Note that Σ_c is not included as an independent parameter since it can be calculated from Equation (27). Also note that the gas and dust components of the disk can each have a different set of values for these parameters.

We first show a fiducial model with stellar and disk parameters listed in Table 2. The disk is assumed to have an inner hole with a sharp edge. Except for the radius of this edge, which is set to 1 AU here, those parameters are set to mimic the transition disk TW Hya as derived by Calvet et al. (2002). The input stellar UV spectrum (including the Lyα line emission) is the observed spectrum of TW Hya (Herczeg et al. 2002, 2004). We did not take into account the possible UV variability. For the initial chemical composition, we assume hydrogen is in H₂, carbon is in gas phase CO, and oxygen not in CO is in water ice, while all the other elements are in atomic form. We assume two types of dust grains, each with a Mathis–Rumpl–Nordsieck (MRN) size distribution (Mathis et al. 1977). The two dust components are assumed to be spatially coexistent in this fiducial model. The larger population has r_min = 1 μm and r_max = 100 μm, with a dust-to-gas mass ratio of 0.01, while the smaller population has r_min = 0.01 μm and r_max = 1 μm, with a dust-to-gas mass ratio of 2 × 10⁻⁵. Larger values for r_max of the big grains has been used in the literature for fitting the disk spectral energy distribution. However, the chemical processes mainly depends on the total available dust grain surface area, which is more sensitive to the assigned overall mass fractions of the small and big grains than the value of r_max. The dust material is assumed to be a 7:3 mixture of "smoothed UV astronomical silicate" and graphite. The optical parameters of the dust are taken from the Web site of Bruce T. Draine⁴ (Draine & Lee 1984; Laor & Draine 1993).

The input density structure, and the gas and dust temperature distribution obtained from radiative transfer and thermal balance calculation are shown in Figure 1. The dominant heating and cooling mechanisms are shown in Figure 2. In the disk upper layer heating is dominated by photoelectric effect, followed by H₂ formation heating in the photodissociation layer, and viscous heating in the deep region. Cooling is dominated by O i and C ii lines in the upper layer, and by accommodation on the dust grains in the lower dense layers. Overall the distribution of the dominant heating and cooling mechanisms is similar to what is shown in Woitke et al. (2009a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The gas phase H₂O distribution is shown in Figure 3. Hot/warm water with temperature ≳200 K is concentrated in a very small region close to the inner edge near the midplane, which can be seen clearly in the inset of this figure. The abundance of water vapor in this region is ≳ 2 × 10⁻⁴, which means essentially all the oxygen not in CO are found in H₂O gas. The total warm water mass is 7.7 × 10²⁶ g, which is equivalent to 560 times the mass of Earth oceans (M_ocean), or ∼0.1 M_Earth. Although the UV flux from the star is very strong (G₀ ∼ 10⁷) at the disk inner edge, the destruction of water by UV radiation is completely quenched only slightly outward, due to the high density in this region. For example, with n_gas = 10¹⁴ cm⁻³ and a dust-to-gas mass ratio of 0.01 (assuming 1 μm for the dust grain radius), the attenuation length of the UV field is of the order of 10⁻⁴ AU, which is too small to be seen in Figure 3. This explains why the warm water is located so close to the inner edge, as envisaged by Cleeves et al. (2011). Similar distributions can also be found in Woitke et al. (2009b), Aresu et al. (2011), Heinzeller et al. (2011), and Meijerink et al. (2012).

Zoom In Zoom Out Reset image size

To understand the small extent of the warm water distribution, it suffices to know that, in a well-shielded region, a high abundance of gas phase water can only be maintained if the dust temperature is higher than ∼150 K, the evaporation temperature of water, otherwise water will condensate onto the dust grains to form ice. The radial temperature gradient close to the inner edge is very steep (see Figure 11) due to scattering and re-emission of radiation toward the disk upper and lower surface, which significantly reduces the amount of energy propagating radially outward. A semi-analytical account of the dust temperature profile in the midplane based on the diffusion approximation of radiative transfer is in Appendix A, where we will see that the dust temperature drops from >300 K at the edge (1 AU) to 100–150 K at r = 1.5 AU, beyond which water can only exist as ice unless some nonthermal desorption mechanisms come into play.

Besides the warm water close to the disk inner edge, there is also cold (≲50 K) gas phase water in the higher layers throughout the disk (main panel of Figure 3). The dust temperature at places where this cold water resides in is well below the evaporation temperature of water, and here the gas phase water is formed from photodesorption of ice (Dominik et al. 2005). The total mass of this diffuse cold water vapor is ∼0.006 M_ocean, but this value depends on the disk size. The overall water vapor mass budget as a function of gas temperature can be seen in the top panel of Figure 4, from which it is clear that a small amount of cold (<80 K) water is associated with reduced but nonzero UV field (indicated by the color scale of the histograms). We note that the total mass of cold water vapor from our model is close to what was derived by Hogerheijde et al. (2011). A small amount of water vapor can also be produced through the radiative association reaction $\mathrm{H + OH \rightarrow H_2O + {\rm h}\nu }$ in the partially photodissociated layer in the outer disk (see also Kamp et al. 2013). Assuming the abundance of water vapor in the outer cold region is determined mainly by photodesorption, photodissociation, and adsorption, and assuming the dust grain is fully covered by water ice, we have

$$\begin{equation} X[\mathrm{H_2O}] = \frac{F_{\rm UV}\sigma _{\rm d} Y \eta _{\rm n}}{F_{\rm UV} \sigma ^{\prime } + n_\mathrm{H} \eta _{\rm n}\sigma _{\rm d}v_{\rm T}}, \end{equation} \tag{ 28 }$$

where F_UV is the UV flux, σ_d is the dust grain cross section, Y is the photodesorption yield, η_n is the dust-to-gas number ratio, σ' is the water photodissociation cross section, and v_T is the thermal speed. As noted by Dominik et al. (2005) and Bergin et al. (2010), when adsorption is unimportant, the abundance of water vapor is independent of the UV flux:

$$\begin{eqnarray} X[\mathrm{H_2O}] &\simeq & \frac{\sigma _{\rm d} Y \eta _{\rm n}}{\sigma ^{\prime }} \nonumber \\ &=& 2.6\,{\times }\,10^{-9}\left(\frac{r_{\rm d}}{0.1\;{\rm \mu m}}\right)^2\left(\frac{\eta _{\rm n}}{10^{-13}}\cdot \frac{Y}{10^{-3}}\right). \quad \end{eqnarray} \tag{ 29 }$$

Deeper into the disk, the density becomes higher and the UV flux becomes weaker, and Equation (29) will overestimate the water vapor abundance with respect to Equation (28).

Also seen in Figure 4 is the existence of a small amount of hot water vapor (≳300 K), which is similar to the result of Woitke et al. (2009b). Except for the region close to the inner wall, where the dust temperature is higher than 300 K, the hot water vapor mainly exists in the upper layer of the disk with r ≲ 30 AU, where the density is low enough that the cooling by accommodation on dust grains is ineffective. The abundance of hot water in such region is ∼10⁻¹⁰, determined by the balance between photodissociation and warm neutral chemistry (see Woitke et al. 2009b; see also Appendix B).

Zoom In Zoom Out Reset image size

In this section, we have seen that the water budget is controlled by a few mechanisms: dissociation by UV photons, adsorption onto dust grains, and the shielding by dust (and possibly self-shielding). Warm water will be preserved in the gas phase if there is enough shielding while keeping the dust temperature higher than the water condensation point. In the next sections, we will discuss their roles in more detail.

We have assumed a razor sharp inner edge in the fiducial model, which gives a well-confined distribution of warm water close to the inner edge. We also run a test model with a "softer" inner edge, namely, we include an exponential taper so that the surface density profile becomes

$$\begin{equation} \Sigma ^{\prime }(r) = \Sigma (r) e^{(r-r_0)/r_{\rm s}},\quad {\rm if }\; r<r_0, \end{equation} \tag{ 30 }$$

where Σ(r) is defined in Equation (26), and r₀ and r_s are the taper parameters. In the test model we let r_in = 0.5 AU, r₀ = 2 AU, and r_s = 0.2 AU, while other parameters are the same as the fiducial model. As shown in the top panel of Figure 5, warm water in this test model is confined in a small region around 1 AU instead of being close to the edge at 0.5 AU simply because the dust shielding only becomes important at ∼1 AU due to the reduced density near the inner edge. Further out from ∼1.5 AU the dust temperature becomes low enough for water to condense out. In this test case, the total mass of warm water is very small, only ∼0.8 M_ocean. There are two reasons for this: the tapered region has much smaller density (hence less mass) than the fiducial one, and the dust shielding in the vertical direction is also reduced (hence smaller volume for water to reside in). The column density of water vapor has a peak value of a few times 10¹⁹ cm⁻² at r ≃ 1 AU, with gas temperature ≳400 K, which are similar to the fitting results of Salyk et al. (2011). For our model to be closer to reality, we will need to solve the disk physical structure self-consistently in a way similar to, e.g., Nomura (2002) or Woitke et al. (2009a), and this will be part of our future study.

The size of the inner hole of a protoplanetary disk may evolve as a result of material exhaustion due to photoevaporation (Gorti & Hollenbach 2009; Owen et al. 2010), or accretion onto the central star or onto the forming planets. Since the warm water is concentrated close to the inner wall as seen in the previous section, we expect the amount of warm water vapor will also evolve as the inner hole expands.

Zoom In Zoom Out Reset image size

We run a set of models with all the parameters except r_in, taking the same value as in the fiducial model. The mass of water as a function of r_in is shown in the top panel of Figure 6. As expected, the warm water mass generally decreases as r_in increases. For r_in > 3.5 AU, the amount of warm water vapor becomes very small because with a central star with bolometric luminosity of 0.25 L_☉, the temperature of the disk wall at 3.5 AU is ∼170 K (taking into account the re-absorption of the radiation emitted by the wall), only marginally higher than the water condensation temperature, hence warm water vapor can only exist in a thin skin of the inner wall.

The bottom panel of Figure 6 shows the vertical column density of water for r_in from 1 to 4 AU. The distribution of water vapor in the inner disk with column densities 10¹⁹–10²² cm⁻² is consistent with the observations of Salyk et al. (2011), though the column densities they derived are mostly in a lower 10¹⁸–10¹⁹ cm⁻² range. The water abundance is ultimately limited from above by the total amount of oxygen available. The narrow rings of water vapor with width ∼0.2 AU resembles what was found by Zhang et al. (2013), though for TW Hya, they apparently found a much higher column density (∼10²² cm⁻², see their Figure 4) compared with 10¹⁹ cm⁻² for the r_in = 4 AU case in Figure 6. One simple way to get a higher column density from our model is to assume a higher surface density at the inner edge, though in our present model the surface density at the edge is already high (∼600 g cm⁻²). One could in principle get a higher oxygen abundance (hence higher water column) relative to hydrogen if the latter has been photoevaporated, without letting the surface density too high. Another possibility to increase the water abundance is related to the details of water adsorption and evaporation (onto and from the dust grains). Since the dust temperature at 4 AU is low (∼160 K), water starts to condense out. The exact temperature for this to happen depends on the dust properties. We have assumed a water desorption energy of 5700 K (Fraser et al. 2001). Lowering this value will release more water into the gas phase. For example, as seen in the gray curve in the bottom panel of Figure 6, reducing it to 5000 K in a test run gives a peak water vapor column density of 10²¹ cm⁻². The desorption energy cannot be too low (≲4000 K) either, because that will tend to overproduce gaseous water, as seen in the black curve in Figure 6. The adsorption and desorption dynamics of dust grains is a complicated issue, which involves the chemical composition, morphology, and crystalline structure of the ice, which themselves are related to their history of formation. The value of 5700 K is appropriate for high density amorphous water ice. Mixing with CO and CO₂ ice (Tielens et al. 1991; Pontoppidan et al. 2008) can reduce the desorption energy of water due to weaker bonding (Cuppen & Herbst 2007), though this is not supposed to happen close to the inner edge of the disk, where CO should be mainly in gas phase. Finally, we must caution that potential degeneracy in parametric fitting of the molecular abundance distribution may render the above comparison premature.

A more luminous central star can warm up the disk, which tends to have more water molecules released into the gas phase, at the same time it also emits more UV photons that can dissociate water molecules (though note that for young stars the UV photons are mainly produced nonthermally by magnetospheric accretion). To explore how the stellar luminosity affect the water distribution, we take a very simple approach by running a set of models with different effective temperature for the central star, assuming the stellar emission to be blackbody without any excess, while keeping all the other parameters the same as in the fiducial model. The mass of warm water vapor as a function of stellar temperature can be seen in the top panel of Figure 7, which shows that the warm water mass increases almost linearly with stellar temperature. The reason is that in our model, a higher stellar luminosity increases the overall temperature of the entire disk. The higher dissociating photon flux associated with a higher T_eff does not reduce the water content because UV photons are readily shielded by the dust, and get converted into photons with lower frequency and no dissociating capability and diffuse into the disk to warm up the dust and gas. The column density profile for different stellar temperatures can be seen in the bottom panel of Figure 7. Higher stellar temperature leads to a wider profile, while the peak column densities are identical, since it is limited by the surface density. Admittedly, our treatment here is simplistic. In reality, a higher overall temperature would "blow up" the disk, reducing the opacity and letting the UV photons penetrate deeper into the disk to dissociate the water molecules, at the same time more water molecules might be liberated into the gas phase due to higher flux of desorbing UV photons.

To see the effect of assuming the stellar spectrum to be a blackbody, Figure 8 shows the column density profile of water vapor in different temperature ranges with (solid line) or without (dashed line) UV excess in the input stellar spectrum, where the stellar temperature is 4000 K and the UV excess is taken to be the same as in the fiducial model. Although a UV excess slightly reduces the column density of the (300–1000) K component, it also increases the column density of the tenuous hot (T > 1000 K) component, for which UV radiation is the major heating source (through photoelectric effect, H₂ photodissociation, H₂ vibrational excitation, etc.).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the disk evolves, dust grains coagulate and grow to larger sizes, and gradually settle down to the midplane (and finally get assembled into forming planets or accreted into the star), leaving less dust in the bulk, and the dust grains that are left at higher altitudes will preferentially have a smaller size than those close to the midplane (but may still be larger than the ISM dust). Our model does not yet contain a self-consistent prescription for dust growth and settling (see, e.g., Dullemond & Dominik 2004b; Dullemond & Spruit 2005), so we simulate the effect of dust settling in a parameterization manner in two ways. The first is to reduce the overall dust mass of the disk, while keeping the dust-to-gas mass ratio equal the fiducial value over the whole disk; the second is to reduce the scale height of the larger grains, keeping a population of small grains well mixed with gas, while the overall dust-to-gas mass ratio stays the same as in the fiducial model.

The resulting warm water distribution as a function of dust-to-gas mass ratio and of the dust scale height can be seen in Figures 9 and 10. The general trend is that less dust present in the bulk of the disk means less warm water vapor. The zoom-in plots in Figure 9 show that the size of the region containing high-abundance water vapor shrinks as the amount of dust is reduced. The main reason is that when dust is reduced, water is more susceptible to UV dissociation. For example, when the dust-to-gas mass ratio is reduced from 0.01 to 10⁻⁴, the UV field strength increase by more than a factor of 10³ in the region where water would otherwise be formed and preserved in gas phase. The self-shielding of water only starts to work at a certain depth, depending on the density and UV intensity, but at such a depth the dust temperature may already drops to a level lower than the water condensation temperature; see the next section for further discussion on this. The reduction of dust does have a small positive effect on water vapor abundance: more water can be retained in the gas phase (if not photodissociated) due to less available adsorption area.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The self-shielding of water can be calculated with

$$\begin{equation} f_{{\rm sh},\mathrm{H_2O}} = e^{-\sigma _\mathrm{H_2O} N_\mathrm{H_2O}}, \end{equation} \tag{ 31 }$$

where $\sigma _\mathrm{H_2O}$ is the photodissociation cross section of water, which is 1.2  ×  10⁻¹⁷ cm⁻² at the frequency of Lyα (van Dishoeck et al. 2006). The dust shielding factor can be calculated with

$$\begin{equation} f_{{\rm sh},\mathrm{dust}} = e^{-\sigma _\mathrm{d} N_\mathrm{d}}. \end{equation} \tag{ 32 }$$

For a silicate grain with radius of 1 μm, the absorption cross section at λ ∼ 0.1 μm is about the same as the geometric cross section, hence we may have

$$\begin{equation} \frac{\sigma _\mathrm{H_2O} N_\mathrm{H_2O}}{\sigma _\mathrm{d} N_\mathrm{d}} \simeq 40 \left[\frac{X(\mathrm{H_2O})}{2\times 10^{-4}}\right]\left[\frac{0.01}{\eta }\right] \left[\frac{r_{\rm d}}{1 \,\mu {\rm m}}\right], \end{equation} \tag{ 33 }$$

where η is the dust-to-gas mass ratio. Hence the absorption due to water can be dominant over dust if water is present in gas phase at its highest possible abundance, as shown by Bethell & Bergin (2009), which is also noted in Ádámkovics et al. (2014).

However, we have seen in the previous section that when the dust is settled or reduced, the amount of water that is present is also reduced. This is because the abundance of H₂ will be reduced due to the enhanced UV field and the fact that less dust surface area is available for its formation when dust is reduced (H₂ can shield itself but nevertheless it needs dust surface to form), which will limit the formation of H₂O from the $\mathrm{O \mathop\rightarrow\limits^{ {H}_2} OH \mathop\rightarrow\limits^{ {H}_2} H_2O}$ chain (Ádámkovics et al. 2014). Also of importance is that even if H₂ is well shielded, CO may still be dissociated because its self-shielding is not as efficient. The produced C atoms can be an important competitor to the above chain to form water. When the gas becomes well shielded by the dust, H₂ becomes the dominant hydrogen bearer and carbon exists as CO, and H₂O can be formed and kept in the gas phase if the dust temperature is higher than the water condensation temperature. Actually, in a few test runs in which we turn off the self-shielding of H₂O in the radiative transfer or we assume all the oxygen not in CO is in gas phase (instead of ice) water for the initial chemical composition, no significant changes in the resulting water vapor mass can be seen, though some small differences are indeed noticeable (see the next paragraph). Hence we may say that the self-shielding of water might be important but only in regions with rather special settings with regard to density structure, UV intensity, and dust abundance. Appendix B contains an approximate semi-analytical account similar to Bethell & Bergin (2009) for the role of water self-shielding.

The role of water self-shielding can also be checked from the output of our code. As already noted in Section 2.1, our models work in an iterative manner. For iteration 0, the initial run, only dust absorption and scattering are included, while water absorption and atomic hydrogen scattering are not included in the radiative transfer. The radiative transfer outputs the distribution of dust temperature and radiation intensity over the whole disk, which is used in the chemical and thermal calculations. The updated chemical composition, specifically, the distribution of H and H₂O, are used in the radiative transfer in the next iteration. This process goes on until the abundance distribution of major species (such as H, H₂, CO, H₂O) do not vary appreciably, though the randomness inherited from the Monte Carlo radiative transfer makes convergence in the usual sense difficult to achieve. Figure 11 shows changes in the radial water vapor abundance profiles at different vertical height as the iteration goes on. The dust-to-gas mass ratio is set to 10⁻⁴, and all the rest parameters are the same as the fiducial model. A close look at this figure shows that the curves for second and third iterations extend inward toward the central star relative to the curve of the first iteration, due to the self-shielding of water. For the midplane where the gas density is very high (10¹⁴ cm⁻³), the relative shift is rather small and almost unnoticeable in the figure. For z = 0.1 and 0.15 AU the effect is more obvious, and the peak abundance of water vapor is also raised by one order of magnitude, similar to what was found by Bethell & Bergin (2009, supporting online material). Higher in the disk atmosphere, with z = 0.2 AU, where the density drops to about 10¹³ cm⁻³, no relative shift can be noticed among different iterations. As calculated semi-analytically in Appendix B (Figure 14), with a density of 10¹³ cm⁻³, a dust-to-gas mass ratio of 10⁻⁴, and a strong UV field, the depth for water to be shielded is ∼0.5 AU. However, at depth 0.5 AU into the disk, the dust temperature will be just low enough (see Figure 1, or maybe more clearly, Figure 12) for water to condense out, and then the water vapor abundance will be determined by photodesorption.