{'fraction_non_alphanumeric': 0.0394371563877851, 'fraction_numerical': 0.028101355397093245, 'mean_word_length': 4.156118439517261, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Water in Protoplanetary Disks: Deuteration and Turbulent Mixing

We investigate water and deuterated water chemistry in turbulent protoplanetary disks. Chemical rate equations are solved with the diffusion term, mimicking turbulent mixing in vertical direction. Water near the midplane is transported to the disk atmosphere by turbulence and destroyed by photoreactions to produce atomic oxygen, while the atomic oxygen is transported to the midplane and reforms water and/or other molecules. We find that this cycle significantly decreases column densities of water ice at r<30 AU, where dust temperatures are too high to reform water ice effectively. The radial extent of such region depends on the desorption energy of atomic hydrogen. Our model indicates that water ice could be deficient even outside the sublimation radius. Outside this radius, the cycle decreases the D/H ratio of water ice from 2x10^-2, which is set by the collapsing core model, to 10^-4-10^-2 in 10^6 yr, without significantly decreasing the water ice column density. The resultant D/H ratios depend on the strength of mixing and the radial distance from the central star. Our finding suggests that the D/H ratio of cometary water (10^-3-10^-4) could be established (i.e. cometary water could be formed) in the solar nebula, even if the D/H ratio of water ice delivered to the disk was very high (10^-2).


INTRODUCTION
Water is one of the most important molecules in star and planet forming regions from the perspective of the oxygen reservoir and the origin of Earth's ocean. It is well established that water ice forms on grain surfaces and is already abundant (∼10 −4 per hydrogen nucleus) in molecular clouds at visual extinctions A V 3 magnitude (e.g., Whittet 1993). Water has also been detected toward star and planet forming regions, such as prestellar cores (Caselli et al. 2012), protostellar envelopes (e.g., Ceccarelli et al. 2000;Boogert et al. 2008), and Class II protoplanetary disks (PPDs; Terada et al. 2007;Honda et al. 2009;Carr & Najita 2008;Hogerheijde et al. 2011). Although there is no doubt that water exists in all the phases of star formation, the evolution of water from molecular clouds to PPDs remains unclear.
D/H ratio can be a useful tool to probe the evolution of water. Comets are preserved samples of the cold ice-bearing regions in the solar nebula (Brownlee 2003). So far, the HDO/H 2 O abundance ratio has been measured in nine comets; seven from the Oort cloud (∼(4-10) × 10 −4 ; e.g., Villanueva et al. 2009;Bockelée-Morvan et al. 2012), and two from the Kuiper belt ((3-<4) × 10 −4 ; Hartogh et al. 2011;Lis et al. 2013). In all sources, the HDO/H 2 O ratio is higher by about one order of magnitude than the elemental abundance of deuterium in the interstellar medium (∼1.5 × 10 −5 ; Linsky 2003). On the other hand, to the best of our knowledge, there has been no detection of HDO toward PPDs in the Class II phase except for . Although they claimed a 4σ detection of HDO vapor toward DM Tau, it has been questioned by Guilloteau et al. (2006).
dimensions. They showed that the majority of water is delivered from the core to the disk as ice. Although they did not include isotopes, their results indicate that D/H ratio of water delivered to the disk is determined in the parent cloud core. In the hot inner regions of disks, where gas temperature is higher than 500 K, isotopic exchange between HDO and H 2 occurs, and the HDO/H 2 O ratio can decrease to ∼10 −5 (Lécluse & Robert 1994). Such reprocessed water could be mixed with unprocessed water via radial mixing and accretion (Drouart et al. 1999;Hersant et al. 2001;Yang et al. 2013). The resultant HDO/H 2 O ratio depends on the thermal history of the disk and efficiency of radial mixing, i.e. how much material diffuse from the inner hot region (>500 K) to outer radii. If we assume the disk formation model of Visser et al. (2011), the HDO/H 2 O ratio decreases by a factor of <2 on average from the initial cloud value, since only ∼30 % of the disk mass experiences high temperatures of several 100 K until envelope accretion onto the protostellar system halts.
The primary motivation of this work is to investigate the impact of turbulent mixing in the vertical direction on the HDO/H 2 O ratio in the Class II PPDs. Deuterium fractionation in Class II PPDs has been studied numerically by several authors without mixing (e.g., Aikawa & Herbst 2001;Aikawa et al. 2002;Willacy 2007;Willacy & Woods 2009). Willacy & Woods (2009) investigated deuterium chemistry in inner disks (≤ 30 AU) considering the accretion flow. In the Class II phase, the disk surface is irradiated by stellar ultraviolet (UV) photons and X-rays. Then, oxygen is mainly in atomic form at the disk surface, while it is in water near the midplane (e.g., Bergin et al. 2007, and references therein). They found that the HDO/H 2 O ratio at the midplane retains their initial value (∼ 2 × 10 −2 , which is set by their dense cloud core model) for 10 6 yr, probably because water is not destroyed efficiently in the regions where UV photons and X-rays are attenuated. However, turbulent mixing could change the situation drastically. If turbulence exists, water would be transported to the disk atmosphere and destroyed by photoreactions, while atomic oxygen would be transported to the disk midplane and reforms water and/or other molecules. Such destruction and reformation processes are of potential importance for the HDO/H 2 O ratio in disks; the HDO/H 2 O ratio could be significantly changed from that in molecular clouds/cores. The plan of this paper is as follows. We describe our physical model and chemical model with vertical mixing in Section 2 and 3, respectively. In Section 4, we present how the vertical mixing affects water and deuterated water chemistry. In Section 5, we discuss the effect of initial HDO/H 2 O ratio, radial mass transport, and dead zone on our results presented in Section 4. We also discuss the uncertainties of grain surface chemistry and the implication of our model results for the D/H ratio of cometary water. We summarize our conclusions in Section 6.

DISK STRUCTURE
We adopt the disk model of Nomura & Millar (2005) with the addition of X-ray heating by Nomura et al. (2007). Similar models have been used to study the chemical evolution without deuterium in protoplanetary disks (Walsh et al. 2010(Walsh et al. , 2012Heinzeller et al. 2011).
Here we briefly describe the disk model. More details can be found in the original papers (Nomura & Millar 2005;Nomura et al. 2007).
The model describes the structure of a steady, axisymmetric Keplerian disk surrounding a typical T Tauri star with the mass M * = 0.5 M , radius R * = 2 R , and effective temperature T * = 4000 K. The α-disk model (Shakura & Sunyaev 1973) is adopted to determine the radial distribution of the surface density, assuming the viscous parameter α = 10 −2 , and accretion rateṀ = 10 −8 M yr −1 . The gas temperature, dust temperature and density distributions of the disk are calculated self-consistently (see Section 2 of Nomura & Millar 2005, for details). The stellar UV luminosity is set to be ∼10 31 erg s −1 in the range of 6-13.6 eV, while the X-ray luminosity is set to be ∼10 30 erg s −1 in the range of 0.1-10 keV based on the spectrum observations of TW Hydra (e.g., Herczeg et al. 2002;Kastner et al. 2002). In the model, the dust-to-gas mass ratio is assumed to be 10 −2 with the dust-size distribution model for dense clouds (Weingartner & Draine 2001). The resultant (a) gas temperature, (b) dust temperature, (c) number density of gases, and (d) wavelength-integrated FUV flux normalized by Draine field (Draine 1978) are shown in Figure 1.

CHEMICAL MODEL
We use a two-phase (gas and grain surface) model to compute disk chemistry (Hasegawa et al. 1992). Our chemical reaction network is based on Aikawa et al. (2012), who investigated deuterium fractionation in the star formation from a prestellar core to a protostellar core. They extended the network of Garrod & Herbst (2006) to include multi deuterated species and isotope exchange reactions from Millar et al. (1989) and Roberts et al. (2004). More detailed explanation can be found in Aikawa et al. (2012). To reduce the computational time, we exclude the species with chlorine and phosphorus, molecules with more than four carbon atoms, and relevant reactions. We also make several modifications to make the network model applicable to disk chemistry, which are described below.
In the regions where the gas temperature exceeds 100 K, we use the reaction set based on Harada et al. (2010) instead of those of Garrod & Herbst (2006) for the neutral-neutral and neutral-ion reactions. The network of Harada et al. (2010) was designed for modeling hightemperature gas phase chemistry (100-800 K). This switching method of the reaction sets was also used in Furuya et al. (2012) to simulate the chemical evolution in collapsing molecular cloud cores. For high temperature reactions of water and its isotopologues, we adopt the rates in Talukdar et al. (1996), Bergin et al. (1999), and the UMIST database (Woodall et al. 2007). They are summarized in Table 2. We include some three-body association reactions and the reverse reactions (i.e. collisional dissociation) based on Willacy et al. (1998) and Furuya et al. (2012). We also modify or include the following reactions, which are important for disk chemistry: (i) X-ray chemistry (Section 3.1 and 3.2), (ii) photoreactions by stellar and interstellar UV (Section 3.2), (iii) chemistry of vibrationally excited H 2 and its isotopologues (Section 3.3). The treatment of gas-grain interactions and grain surface reactions is described in Section 3.4. In total, our reaction network consists of 726 gaseous species, 303 grain-surface species, dust grains with three different charge states (±1, 0), and 34,997 (38,162) reactions for gas temperatures ≤100 K (>100 K).

Ionization Rate
Interstellar cosmic-rays and X-rays from a central star are the main ionization sources of protoplanetary disks. The cosmic-ray ionization rate of hydrogen molecules is expressed as where Σ 1 (r, z) and Σ 2 (r, z) are the vertical column density of gas measured from the upper and lower boundary of the disk to the height z, respectively (e.g., Semenov et al. 2004). Unattenuated cosmic-ray ionization rate, ξ 0 CR , is set to be 5×10 −17 s −1 (Dalgarno 2006). The attenuation length of cosmic-ray is 96 g cm −2 (Umebayashi & Nakano 1981), which is larger than our disk column density only at r 1 AU. We also take into account the ionization by short-lived radio active isotopes with the rate of 1 × 10 −18 s −1 (Umebayashi & Nakano 2009).
The primary ionization rate of species i by X-rays is given by where F X is the photon number flux from the central star. The photoionization cross sections, σ pi,i , are taken from Verner & Yakovlev (1995). We consider the primary ionization of atoms, singly ions, and diatomic molecules. Considering the Auger effect, it is assumed that ionization of atoms leads to doubly charged ions, while ionization of diatomic molecules produces a pair of singly charged ions, following Maloney et al. (1996). The cross sections of the molecules are calculated by adding the atomic cross sections. We use the same cross sections for deuterated species as normal isotopes. Electron recombination and charge transfer reactions for doubly ionized species are taken from Stäuber et al. (2005). We neglect doubly ionized states of Na and Mg, since the elemental abundances of Na (2.25 × 10 −9 ) and Mg (1.09 × 10 −10 ) are low in our model. We also neglect doubly ionized state of Si, since its abundance is much lower than that of Si + in PPDs (e.g., Ádámkovics et al. 2011).
Part of the kinetic energy of photoelectrons produced by primary ionization is lost by ionizations of H, H 2 , and He. These secondary ionizations are far more important for the ionization degree of gas than the primary ionization. Assuming that the kinetic energy of Auger electrons is equal to the photoionisation threshold energy (Voit 1991), the overall secondary ionization rate per hydrogen nucleus is where σ eff (= i x cosmic i σ pi,i ) is the effective photoionization cross section per hydrogen nucleus evaluated assuming a cosmic abundance for elements (x cosmic i ). Dalgarno et al. (1999) derived the mean energy required per secondary ionization of species i, W i , for the H-H 2 -He mixtures as a function of the ionization degree and the abundance ratio of H and H 2 . We use Equations (9), (10), (13) and (15) of Dalgarno et al. (1999) with the parameter applicable to a 1 keV electron to calculate W in our models (e.g., Meijerink & Spaans 2005). Although secondary ionizations of other species have no significance for the ionization degree of gas, it could be important for chemistry (Maloney et al. 1996). The secondary ionization rate of other species i per hydrogen nucleus is estimated by where x i represent the abundance of species i relative to hydrogen nuclei. The electron impact cross sections, σ ei , are taken from the NIST database (http://www.nist.gov/pml/data/ionization/) and Lennon et al. (1988). We include secondary ionization of CO, singly charged atomic ions, and atoms except for H and He in our network using Equation (4). In our disk models, the secondary ionization by X-rays is the dominant ionization source except for the midplane (e.g., N H 10 24 cm −2 at 5 AU, where N H is the vertical column density of hydrogen nucleus from the disk surface.), where cosmic ray ionization is dominant.

Photochemistry
Protoplanetary disks are irradiated by UV photons from interstellar radiation and central stars. While the former has a continuous spectral energy distribution, the latter is the combination of continuum and Lyα radiation. Furthermore, cosmic-rays and X-rays induce UV photons through energetic photoelectrons. The UV spectrum at each point in the disk is given as a combination of them. Because of the different origins and wavelength dependence, photochemistry by these photons is considered separately. In this subsection, rate coefficients (s −1 ) are denoted by k, while photorates (cm −3 s −1 ) are denoted by R.

UV Continuum
Spectral energy distribution of UV continuum is given at each point in the disk (Section 2). To calculate photorates, we convolve the UV spectrum with absorption cross sections. Referring to van Dishoeck (1988), the rate coefficients of photodissociation and photoionization in the gas phase (s −1 ) are calculated as where J(λ) is the the mean intensity of the radiation field measured in photons cm −2 s −1 A −1 str −1 at the wavelength of λ. The cross sections, σ(λ), are given by van Dishoeck (1988), and updated by Jansen et al. (1995a,b) and van Dishoeck et al. (2006). The data is downloadable from http://www.strw.leidenuniv.nl/˜ewine/photo/. For deuterated species, we use the same cross sections as normal species. For species which are not in the database, we use the cross sections of a similar type of species (e.g., van Zadelhoff et al. 2003).
H 2 , HD, and CO are dissociated by absorption lines. Because of their high abundances, self-shielding and mutual shielding of HD and CO by H 2 are important. The selfshielding factor of H 2 is given as a function of H 2 column density referring to Equation (37) of Draine & Bertoldi (1998). Since the self-shielding behavior of HD is almost identical to that of H 2 (Wolcott-Green & Haaiman 2011), we also use Equation (37) of Draine & Bertoldi (1998) for self-shielding factor of HD by replacing H 2 column density with that of HD. The shielding factor of HD by H 2 is given by Equation (12) of Wolcott-Green & Haaiman (2011). For CO, the self-shielding and mutual shielding factors are calculated following Lee et al. (1996).
Then, the photodissociation rate coefficient of CO, for example, is calculated as follows: where k 0 phg,r and k 0 phg,z are non-shielding photodissociation rates of CO, which are calculated at each point using Equation (5), by the UV radiation in the radial and vertical directions, respectively. The shielding factor, f sh , is a function of the column densities of CO and H 2 , N CO and N H 2 , along the ray. In the vertical direction, we obtain the shielding factor at each vertical position z by integrating the number density of CO and H 2 at >z, where n i represents the number density of species i. In the radial direction, on the other hand, molecular column densities are assumed to be the overall column density along the ray multiplied by the local molecular abundances: where n H represent the number density of hydrogen nuclei. Equation (8) could overestimate the column densities (i.e. shielding factors), since the radiation field decays and the abundances of CO and H 2 would increase along the ray.
Grain surface species are also photodissociated and photodesorbed by UV photons (e.g., Gerakines et al. 1996;Westley et al. 1995). The photoabsorption cross sections of H 2 O ice and CO 2 ice are taken from Mason et al. (2006). It is worth noting that the cross sections of ices are different from those in the gas phase (see Figure 1 of Andersson & van Dishoeck 2008, for water). For the other species, we use the same cross section for corresponding gaseous species, because data is not available in the literature. We confirmed that our results are not sensitive to this assumption. We performed a calulation using the ten times smaller cross sections for ice molecules except for water and CO 2 ices, and found that the results are not significantly different from those in Section 4. The photodissociation rates on grain surfaces (cm −3 s −1 ) are calculated as R phs,i (r, z) = πa 2 n d θ i 4πJ(λ, r, z)P i (λ)dλ, where n d , n ice , θ i , and, N site are the number density of the grain, the total number density of grain surface species, the coverage of species i on a grain, and the number of adsorption sites on a grain (N site ∼ 10 6 ), respectively. The grain radius, a, is set to be 0.1 µm in our chemical model for simplicity. The ratio of the geometrical and absorption cross sections of dust grains for FUV photons is assumed to be unity (e.g., Tielens 2005). P i (λ) and σ site are the probability for an incident photon with the wavelength of λ on a site to be absorbed by an adsorbed species and the site size (4πa 2 /N site ), respectively.
Every photoabsorption does not necessarily result in photodissociation. Andersson et al. (2006) and Andersson & van Dishoeck (2008) carried out molecular dynamics simulations to study photodissociation of water ice in the top six monolayers. They found that photolytic effect is more significant in the uppermost few monolayers; although photons can penetrate into deeper layers and dissociate embedded molecules, the photoproducts recombine with a higher probability than in uppermost layers. In this study, we assume that the upper most one monolayer can be dissociated as the outcome of photoabsorption, and that the photoproducts immediately recombine in deeper layers. For considering surface roughness, parameter N p in Equation (9) is set to be two monolayers for all species. We discuss the effect of this assumptions on our results in Section 5.2.
In this study, we adopt a two-phase model, in which a layered structure of ice mantles is not considered and the bulk of ices is chemically active. More sophisticated ice mantle modeling, which distinguish a chemically active surface layer from an inactive inert bulk, has been developed (Hasegawa & Herbst 1993;Taquet et al. 2013b). Since the calculation of the rate equations combined with turbulent mixing is time consuming, we postpone the inclusion of such layered ice structure to future works. The combination of the two-phase model and photodissociation of HDO ice could, however, introduce artificial decrease of the D/H ratio of water ice as follows. Photodissociation of HDO ice has two branches, producing OH ice or OD ice. While OH ice is able to cycle back to HDO ice via the grain surface reactions, e.g., OH ice + D ice , it is also converted to H 2 O ice. If photodissociation of ices is efficient only in the uppermost several layers, this cycle could decrease the D/H ratio of water ice on the surface, but not the ratio in the bulk water ice. In the two phase model, on the other hand, this cycle could decrease the D/H ratio in the bulk water ice; we do not distinguish between D/H ratios in the surface layers and the deeper layers, although we calculate the photodissociation rate of ice assuming that only the uppermost one layer can be dissociated. To keep the consistency in our model, we switch off the branch to produce OH ice in photodissociation of HDO ice.
The photodesorption rates are calculated as (e.g., Visser et al. 2011) where F FUV and Y i are the FUV photon number flux integrated in the range of 912-2000Å and the photodesorption yield per incident FUV photon, respectively. Some experiments into the photodesorption of UV-irradiated ices have been conducted (e.g., Westley et al. 1995;Oberg et al. 2007). Recently,Öberg et al. (2009a,b) measured the yield of CO 2 and H 2 O ices by photons with the wavelength of 1180-1770Å. Fayolle et al. (2011Fayolle et al. ( , 2013) measured wavelength-dependent photodesorption yields of CO, O 2 , and N 2 ices, and calculated the yields adequate for astrophysical environments. For the above species, we use the yields obtained by the experimental works. The yield for the other species is assumed to be 10 −3 . We confirmed that the results presented in this paper does not change, if we use the yield of 10 −4 for species without laboratory data.

Lyα Photons
Lyα radiation could dominate the UV field in protoplanetary disks (Herczeg et al. 2004). We include photoreactions by Lyα photons for a limited number of gaseous species, for which photodissociation cross section at Lyα wavelength (∼1216Å) is available in van Dishoeck et al. (2006), and corresponding species in the solid phase. The cross sections for H 2 O ice and CO 2 ice are taken from Mason et al. (2006), while we use the cross section of the gaseous species for the other ice-mantle species. Considering that the photodissociation of ice-mantle species causes photodesorption (Andersson & van Dishoeck 2008), the photodesorption yields per incident Lyα photon, Y Lyα , are estimated by where σ i (λ Lyα ) andσ i are the photodissociation cross section at Lyα wavelength and the averaged cross section in the wavelength of 1180-1770Å, respectively.
Recently, Bethell & Bergin (2011) showed that the Lyα/FUV-continuum photon density ratio is smaller in the regions near the photodissociation layers for hydrogen (N H 5 × 10 20 cm −2 at 1 AU if the dust is well mixed with gas), but larger in the deeper regions than the intrinsic ratio in the incident FUV radiation from the central star, due to resonant scattering of Lyα by atomic hydrogen. Since the intensity of Lyα radiation is calculated without considering resonant scattering in the physical model used here, it would be overestimated in the upper layers of the disk, and underestimated near the midplane.

Cosmic-ray and X-ray Induced UV Photons
UV photons are generated by the decay of electronically excited hydrogen molecules and atoms produced by collisions with energetic photoelectrons following cosmic-ray ionization. In dense regions where the stellar and interstellar UV photons are heavily attenuated, cosmicray induced photons dominate the UV fields (Prasad & Tarafdar 1983). X-ray ionization also induces UV photons in a similar manner. Our model includes photochemistry by cosmicray and X-ray induced UV photons both in the gas phase and on grain surfaces. The rates are described by the local competition of photon absorption by gaseous species and dust grains. Referring to Gredel et al. (1989), the rate coefficients of photodissociation and photoionization of species i in the gas phase (s −1 ) are expressed as where ψ(λ), and x d are the number of photons with the wavelength of λ produced per ionization, and the relative abundance of dust grains to hydrogen nucleus (∼ 2 × 10 −12 ), respectively. Again, the ratio of the geometrical and absorption cross sections of dust grains for FUV photons is assumed to be unity (e.g., Tielens 2005). For simplicity, we ignore photons produced by the de-excitation of atomic hydrogen and assume x H 2 = 0.5 in Equation (14). It is an acceptable assumption, since the regions dominated by atomic hydrogen do not contain many molecules (Aikawa & Herbst 1999). We also ignore the gas opacity in Equation (15). If all oxygen is in gaseous water, the latter assumption overestimates the rate coefficients by a factor of about two at maximum, since However, it is also acceptable, since induced UV photons are more important than stellar UV photons only near the midplane, and dust temperatures near the midplane in our disk models are lower than the sublimation temperature of water ice except for r 1 AU.
The rates (cm −3 s −1 ) of photodissociation on grain surfaces and photodesorption are calculated by Figure 2 shows the timescales of photoreactions of water ice and water vapor as a function of vertical column density in our disk model at r = 10 AU. Here, the photoreaction timescales of water ice are defined as the timescales in which one monolayer of pure water ice is photodissociated (photodesorbed), n d N site /R(θ H 2 O ice = 1), while those of water vapor are defined as an inverse of the rate coefficients, k −1 . It is clear that the photoreactions by Lyα photons are most efficient in the disk atmosphere, while the photoreactions by X-ray and cosmic-ray induced UV photons are the most efficient at N H 2 × 10 22 cm −2 at 10 AU.

Photoreaction Timescales in the Disk
It should be noted, however, that these timescales are larger than the typical age of T Tauri stars (∼10 6 yr) at N H 10 23 cm −2 .

Vibrationally Excited Hydrogen Molecule
Hydrogen molecules can be electronically excited by absorbing FUV photons. Subsequent fluorescence leads to dissociation or to the population of various rotation-vibration levels in the ground electronic state. The branching ratio of the former and latter are 10 % and 90 %, respectively (Black & Dalgarno 1976). Following Tielens & Hollenbach (1985), we consider only two states of H 2 , the ground vibrational state and a single vibrationally excited state (H * 2 ). They are treated as different species in our chemical model. The effective quantum number of this pseudolevel is v = 6, and the effective energy is 2.6 eV (∼3 × 10 4 K, London 1978). The FUV pumping rate is set to be 9 times the photodissociation rate of H 2 . H * 2 can be destroyed by photodissociation, radiative decay, collisional de-excitation, or collisional dissociation. These rates are calculated following Tielens & Hollenbach (1985). As an inverse reaction of the collisional de-excitation, we consider collisional excitation following Woitke et al. (2009a).
H * 2 reacts with other species using internal energy to overcome activation barriers. It should be noted, however, that whether the internal energy of H * 2 is effectively used or not is quite specific for each reaction and difficult to predict as discussed by Agúndez et al. (2010). Therefore, it would be problematic to simply assume that all the internal energy of H * 2 is used to overcome activation barriers in any reactions. To be conservative, we include a limiting number of reactions between H * 2 and other species (He + , C + , O, OH, and CN), rate coefficients of which are experimentally or/and theoretically determined (see Table 1 of Agúndez et al. 2010). Vibrationally excited HD and D 2 are treated in the same way as H * 2 .

Gas Grain Interactions and Surface Reactions
Calculations of gas-grain interactions and grain-surface reactions except for H 2 formation are performed in a similar way to Garrod & Herbst (2006) and Furuya et al. (2012). A brief summary is provided below.
The sticking probability of neutral species onto dust grains is assumed to be unity except for atomic hydrogen and deuterium; their sticking probabilities are calculated as a function of gas and dust temperatures following Hollenbach & Mckee (1979). Interactions between ions and dust grains are calculated in the same way as Furuya et al. (2012). We adopt the same desorption energies (E des ) of atoms and molecules as Garrod & Herbst (2006), in which grain-surfaces are assumed to be covered with non-porous water ice, and thus only physisorption sites are considered, unless stated otherwise. For deuterated species, we use the same desorption energies as normal species except for atomic deuterium, whose adsorption energy is set 21 K higher than that of atomic hydrogen, following Caselli et al. (2002). E des of selected species are listed in Table 3. In addition to the thermal desorption and photodesorption (see Section 3.2), we consider sublimation via stochastic heating by cosmic-ray (Hasegawa & Herbst 1993) and X-ray, and the chemical desorption. We assume that roughly 1% of species formed by surface reactions are desorbed following Garrod et al. (2007).
Grain surface reactions are assumed to occur by the Langmuir-Hinshelwood (LH) mechanism between physisorbed species; adsorbed species diffuse by thermal hopping and react with each other when they meet (e.g., Hasegawa et al. 1992). The energy barrier against diffusion is set to be a half of the desorption energy. If surface reactions have activation energy barriers, they are overcome thermally, or via quantum tunneling, whichever is faster (Garrod & Herbst 2006). We use the modified rate method (Caselli et al. 1998) for surface reactions which include atomic hydrogen or deuterium as reactants.
In our disk models, H 2 formation by association of physisorbed H atoms does not efficiently occur by the LH mechanism, since the dust temperatures are 20 K except for the midplane at r 250 AU (e.g., Hollenbach & Salpeter 1971). Cazaux & Tielens (2004, 2010 and Cazaux et al. (2008) studied formation of H 2 and its isotopologues on bare grain surfaces, considering both physisorbed and chemisorbed sites and both the LH and Eley-Rideal mechanisms. They showed that H 2 and HD can be formed with the efficiency of several tens of percent until quite high temperatures (several hundreds K). Since the abundance of water in disk atmospheres is highly dependent on the abundance of H 2 (e.g., Glassgold et al. 2009), H 2 formation on bare grain surfaces should be considered. The precise efficiency is, however, highly dependent on the width and height of the barrier between the adsorption sites, which is not well constrained; we adopt a constant efficiency of 0.2 for H 2 and HD independent of dust temperature. We ignore D 2 formation on bare grain surfaces, since it is not as efficient as that of H 2 and HD (Cazaux et al. 2008). Considering the accretion rate of a pair of atomic hydrogen onto bare grain surfaces, the formation rate of H 2 is expressed as where ε H 2 is the formation efficiency. The coverage of ice mantle species, θ ice , is given by substituting n ice into n i in Equation (10). The rate is set to zero when dust grains are covered by more than one monolayer of ice (i.e. n ice ≥ n d N site ).

Initial Abundances for Disk Chemistry
We adopt the so-called low metal values as the elemental abundances (see Table 1 of Aikawa & Herbst 2001). The elemental abundance of deuterium is set to be 1.5 × 10 −5 (Linsky 2003). It is assumed that hydrogen and deuterium are initially in H 2 and HD, respectively. The heavy elements are assumed to be initially in atomic or ionic form, corresponding to their ionization energy. We integrate the rate equation, using our chemical network model, in the collapsing core model of Aikawa et al. (2012), in which fluid parcels are traced from the prestellar core to the protostellar core of age 9.3 × 10 4 yr. Because the abundances are mostly constant at r 100 AU in this protostellar core model, we adopt the molecular abundances at r = 60 AU as the initial abundance of our disk model (Table 4).

Turbulent Mixing
We compute molecular evolution in the protoplanetary disk with vertical mixing, assuming that the origin of the disk turbulence is the magnetorotational instability (MRI; Balbus & Hawley 1991). We consider the inhomogeneous one-dimensional diffusion equations at a specific radius (e.g., Xie et al. 1995;Willacy et al. 2006), where P i , and L i represent the production rate and the loss rate, respectively, of species i. The second term in the left-hand side of Equation (20) describes diffusion in turbulent disks. The vertical diffusion coefficient for gas and very small dust grains (i.e well-coupled to gas) is assumed to be where Ω is the Keplerian orbital frequency (Fromang & Papaloizou 2006;Okuzumi & Hirose 2011). The vertical velocity dispersion, δv 2 z , is assumed to be where c s is the local sound speed. The Schmidt number (ratio of gas to dust diffusivity) is where τ s is the stopping time of dust grains (Youdin & Lithwick 2007). We adopt Sc = 1, since the term Ωτ s is much less than unity in the midplane of our disk model (e.g., ∼10 −6 at r = 5 AU). In other words, gaseous species and dust grains (grain surface species) have the same diffusivity. It is worth noting that in our approach, only the mean composition of ice mantle of grains is obtained at each spatial point in the disk. In reality the ice composition could vary among grains in the same spatial grid at a given time, since the motion of grain is random in turbulent gases, and since grains have different thermal and irradiation history (Ciesla & Sandford 2012).
The activity of MRI depends on the ionization degree. MRI can be stabilized near the midplane of disks, where the ionization degree is low and magnetic field is decoupled to plasma. It is called dead zone (e.g., Gammie 1996;Sano et al. 2000). Three-dimensional isothermal MHD simulations in the shearing box approximation show δv 2 z in the dead zone is smaller than that in the MRI active layer by more than one order of magnitude, depending on the strength of the magnetic field (Okuzumi & Hirose 2011;Gressel et al. 2012). In our fiducial models, however, we assume that α z is constant for simplicity. We run three models in which α z is set to be 0, 10 −3 , and 10 −2 . The effect of dead zone is discussed in Section 5.4.
Equation (20) is integrated for 10 6 yr using implicit finite differencing on a linear grid consist of vertical 60 points at a specific radius. Our code is based on Nautilus code (Hersant et al. 2009), but we improved the treatment of turbulent diffusion to account for the full coupling of mixing with chemistry according to Heinzeller et al. (2011). As boundary conditions, we assume that there is no flux through the midplane and the upper boundary of the disk. We integrate Equation (20) at 35 radial points from r = 1 AU to r = 300 AU. We ignore radial accretion and radial mixing in the present study; we will investigate them in forthcoming papers.

Effect of Vertical Mixing on Oxygen Chemistry
Beyond the snow line (r 1 AU in our models), the bulk of water exists as ice on grain surfaces. In this subsection, we describe how the vertical mixing affects abundances of oxygen reservoirs beyond the snow line in the disk. As an example, we look into the results at r = 10 AU and 30 AU in the model with α z = 10 −2 . In the following, species X on grain surfaces are denoted as X ice . Figure 3 shows physical parameters (top panels), abundances of assorted O-bearing species in the model with α z = 0 (middle panels), and α z = 10 −2 (bottom panels) as functions of vertical column density at r = 10 AU (left panels) and 30 AU (right panels), at t = 10 6 yr. In the model without mixing, the disk can be divided into two layers in terms of the main oxygen reservoir; atomic oxygen in the disk surface and water ice near the midplane. Above the O/H 2 O ice transition, photoreactions on grain surfaces are efficient, and the water ice abundances are low (<10 −6 ). Thermal desorption of water ice is not efficient either at r = 10 AU or r = 30 AU, since dust temperatures are less than the sublimation temperature of water ice (∼150 K) even at the disk surface.
In the model with mixing, water ice is transported to the surface through the O/H 2 O ice transition. Such water ice is destroyed by photoreactions to produce atomic oxygen. On the other hand, atomic oxygen is transported to the deeper layers in the disk, and is mainly converted to O 2 at r = 10 AU, and CO 2 ice and water ice at 30 AU. Then, the water ice abundance near the midplane decreases with time; it is less than the canonical value of ∼10 −4 at 10 6 yr especially at r = 10 AU. At r 40 AU, on the other hand, water ice abundance is ∼10 −4 even in the model with mixing (see Figure 5), since most of atomic oxygen is converted back to water ice. The variation of the major O-bearing species, which are (re)formed from the atomic oxygen, partly depends on dust temperatures near the O/H 2 O ice transition (strictly speaking the height z * ; see Section 4.2), where the conversion of atomic oxygen to other species mainly takes place. In Figure 4 we present main formation paths to (re)form O 2 , CO 2 ice, and water ice beyond the snow line in our models. We note here that the formation reactions of these species include OH as the reactant: Reaction (25) is a gas-phase reaction, while the others are grain-surface reactions. In order for water ice to be efficiently reformed, OH should be mainly formed on grain surfaces, and water ice formation should be more efficient than CO 2 ice formation.
At r = 10 AU, dust temperatures are too high (∼60 K) to effectively form OH on grain surfaces; adsorbed atomic oxygen is likely to be evaporated before it meets other reactive species. Hence, the main formation route of OH from atomic oxygen is the gas phase reactions: radiative association with atomic hydrogen and/or a sequence of ion-neutral reactions, followed by the recombination of H 3 O + with an electron. In the recombination of H 3 O + with an electron, the branch to produce OH accounts for 74 % of the total recombination rate (Jensen et al. 2000). Although the recombination of H 3 O + also has a branch to produce water vapor, it accounts only for 25 % of the total rate. Once OH is formed in the gas phase, it mostly reacts with atomic oxygen and is converted to O 2 before being adsorbed onto dust grains, due to the high abundance of atomic oxygen ( 10 −5 ) around the O/H 2 O ice transition. The ratio of the O 2 formation rate and adsorption rate of OH is expressed as where k 1 , the rate coefficient of OH + O, is 4 × 10 −11 cm 3 s −1 independent of the gas temperature. Then, the conversion to O 2 is dominant. Although the value of k 1 has some uncertainties (2-8×10 −11 cm 3 s −1 ) according to KIDA (http://kida.obs.u-bordeaux1.fr; see also Wakelam et al. 2012;Hincelin et al. 2011), the ratio exceeds unity even if we use the lowest value. We note that conversion of O 2 to water on grain surfaces (O 2ice (Miyauchi et al. 2008) is not efficient due to high dust temperature. Another reason for the high abundance of O 2 is that its destruction rate by UV photons is small compared to that of water ice, since Lyα cross section of O 2 is lower than that of water ice by about two orders of magnitude Mason et al. 2006).
At r 30 AU, lower dust temperatures ( 40 K) allow OH to mainly form on grain surfaces. Water ice is mainly formed by H ice + OH ice at r 30-40 AU. CO and atomic hydrogen are mostly in the gas phase at these radii. Since CO has a higher desorption energy than atomic hydrogen (Table 3), adsorbed CO stays for a longer time on grain surfaces andCO ice + OH ice proceeds faster than H ice + OH ice . CO 2 ice has a smaller photoabsorption cross section than that of water ice, which further ensures the high abundance of CO 2 ice at r ∼ 30 AU. H 2 CO has a higher desorption energy than CO, and mostly exists as ice at r 30 AU. At r 40 AU, H 2 CO ice + OH ice is more efficient than H ice + OH ice and CO ice + OH ice , due to high abundance of H 2 CO ice near the O/H 2 O ice transition ( 10 −8 ). Figure 5 shows two-dimensional distributions of the abundance of water ice, while Figure  6 shows radial distributions of its column density. In the model without mixing, the water ice abundance beyond the snow line stays constant. In the model with mixing, the water ice abundance near the midplane decreases. The tendency is more significant at smaller radius, since reformation of water ice is prohibited by high dust temperatures, and since the timescale of mixing is small (see Section 4.2 for quantitative discussions).
Above the O/H 2 O ice transition, a fraction of oxygen exists as water vapor, and its abundance reaches a maximum near the H/H 2 transition. In the regions with the gas temperature T 300 K, water vapor is formed via the neutral-neutral reaction which has the potential energy barrier of 1740 K. Water vapor is also formed by a sequence of ion-neutral reactions, followed by the recombination of H 3 O + with an electron, which is the dominant formation path in the regions with gas temperature T 300 K.
In the model without mixing, we can see the local bumps of the abundances of atomic hydrogen at N H ∼ 10 22 cm −2 . The positions correspond to where dust grains start to be covered with ice and chemisorption sites become unavailable (i.e. formation rate of hydrogen molecule on grain surfaces decreases). These bumps are smoothed out via mixing. In the model with mixing, peak abundance of water vapor is enhanced in the disk atmosphere, since the H 2 abundance above the H/H 2 transition increases (Reaction 30; see also Heinzeller et al. 2011).

Destruction Timescale of Water Ice
When turbulence exists, water ice near the midplane is destroyed via the combination of vertical mixing and photoreactions. This process significantly changes oxygen chemistry as mentioned above. The destruction timescale of water ice would be determined by the timescale of vertical transport via mixing, which in turn, can be estimated from the column density of water ice divided by the upward flux φ H 2 O ice (z * ): We define z * as the height where ∂φ/∂z reaches a minimum (i.e. where the net upward flux of water ice reaches a maximum). We note that the precise evaluation of the timescale requires the solution of Equation (20), since the flux φ is proportional to the composition gradient. We evaluate τ mix H 2 O ice using our numerical data at t = 10 4 yr in the model with mixing. The upper panel of Figure 7 shows temporal variation of water ice abundance in the model with α z = 10 −2 at the midplane of r = 5 AU and 10 AU. It shows that the water ice abundance at the midplane at t < a few 10 5 yr in the numerical simulations are well reproduced by the exponential decay: It ensures that τ mix H 2 O ice represents the characteristic timescale of destruction of water ice. At t > a few 10 5 yr, the water ice abundance is significantly low compared to the initial value, and destruction and formation almost balance. The lower panel of Figure 7 is similar to the upper panel, but shows the results at r = 30 AU and 50 AU. At these radii, Equation (32) is no more valid, since water ice reformation is not negligible, although Equation (31) should be still valid.
We can estimate the destruction timescale of HDO ice from Equation (31), but with the subscript replaced by HDO ice. We find that the timescales of HDO ice and H 2 O ice are similar, because their chemistry is similar. For example, τ mix H 2 O ice and τ mix HDO ice are 6.3 × 10 4 yr and 5.7 × 10 4 yr, respectively, at r = 10 AU. Then we simply denote them as τ mix in the rest of the paper.
In Figure 8, we present radial variations of τ mix . We fit the numerical data at 2 AU r 100 AU and find that τ mix depends almost linearly on the radial distance from the central star: τ mix (r) ∼ 6 × 10 4 r 10 AU α z 10 −2 −0.8 yr (2 AU r 100 AU).
The timescale is smaller in the inner regions. At r 100 AU, destruction timescale is not determined by τ mix , since photoreactions rather than the upward transport of water ice limits the destruction rate, which is clear from the flat distribution of water ice abundance in the vertical direction ( Figure 5).

Deuteration of Water Ice
We describe how the vertical mixing affects deuterium chemistry, focusing on the largest deuterium reservoir, HD. In our initial abundance for the disk model, which is set by the collapsing core model, ∼30 % of deuterium is in species other than HD and D 2 (Table 4). For example, ∼15 % of deuterium is in water, and its D/H ratio is ∼2 × 10 −2 . In the collapsing core, the species which are mostly formed in a cold era (∼10 K), such as water, are enriched in deuterium through the isotope exchange reactions, such as, H + 3 + HD ⇋ H 2 D + + H 2 + 230 K. The backward reaction is endothermic, and is negligible compared to the forward reaction in the low gas temperature. Then, H + 3 and several other species which are subject to the direct exchange reactions become enriched in deuterium. Chemical reaction network propagates the enrichment to other species, so that a significant amount of deuterium is delivered to the species other than HD and D 2 (e.g., Aikawa et al. 2012). Figure 9 shows vertical profiles of a fraction of elemental deuterium in the form of HD and D 2 (i.e. (x HD + 2x D 2 )/1.5 × 10 −5 ) in the disk models at r = 30 AU. The total fraction in other molecular form than HD and D 2 , mainly ice, are also shown in Figure 9. At N H ∼ 10 20 cm −2 , these fractions are small, since deuterium is mostly in atomic form. In the model without mixing, almost all deuterium is in HD (and D 2 ) at N H = 10 21 -10 22 cm −2 , since the other major deuterium reservoirs, such as HDO ice, are mostly destroyed by photoreactions. In the deeper regions, on the other hand, those reservoirs survive for 10 6 yr, since the timescales of photoreactions are longer than 10 6 yr there (see Figure 2). In the model with mixing, the deuterated ices are transported to the surface to be destroyed by photoreactions, while HD is transported to the midplane. Note that the total elemental abundances in gas and dust, are not changed by the mixing, as long as dust grains are dynamically well-coupled to the gas. Because of the higher gas temperatures in the disk than that in the cold core, the endothermic exchange reactions, such as H 2 D + + H 2 , are not negligible. Therefore, reformation of deuterated ices is not efficient enough to compensate for the destruction. As a result, the HD abundance near the midplane increases with time, while the abundances of the other deuterium reservoirs decrease. For example, ∼20 % of deuterium is in the species other than HD at t = 10 5 yr in the midplane of r = 30 AU, while the fraction decreases to less than 1 % at t = 10 6 yr.
As shown in the bottom panel of Figure 7, the H 2 O ice abundances are almost constant in the midplane at r = 30 AU and 50 AU, since reformation of H 2 O ice (partly) compensates for the destruction. On the other hand, the HDO ice abundances decrease in a similar timescale to τ mix . Then, the HDO ice /H 2 O ice ratio decreases in a similar timescale to τ mix at these radii. In the inner radii, r = 5 AU and 10 AU, reformation of both ices is negligible, due to high dust temperatures. Since the destruction timescales of H 2 O and HDO ices are similar, substantial changes of the HDO ice /H 2 O ice ratio does not occur until their abundances decrease significantly and reach the steady state (t > a few 10 5 yr; upper panel of Figure 7). Therefore, τ mix gives the approximate timescale of the decline of the HDO ice /H 2 O ice ratio at r 30 AU in our models, while it gives the timescale of H 2 O ice decrease at r 30 AU. Figure 10 shows radial variations of the HDO ice /H 2 O ice column density ratio. While the D/H ratio retains its initial value (∼2 × 10 −2 ) for 10 6 yr in the model without mixing, it decreases with time in the model with mixing. Since τ mix is smaller and the gas temperature is higher in inner radii, the resultant D/H ratio of water ice column density is roughly an increasing function of radius. In the model with α z = 10 −3 , the HDO ice /H 2 O ice ratio decreases by more than one order of magnitude at r 3 AU within 10 6 yr. In the model with more efficient mixing, α z = 10 −2 , such a significant decline of the HDO ice /H 2 O ice ratio can be seen inside r = 2 AU and 30 AU, within 10 5 yr and 10 6 yr, respectively. It should be noted, however, that in these regions with the low HDO ice /H 2 O ice ratio, water ice abundance and its column density are also lower than those in the model without mixing by about two orders of magnitude, and water ice is no longer the major oxygen reservoir.
In summary, at r < 30 AU, warm gas temperatures lead to low HDO ice /H 2 O ice ratio, while warm dust temperatures hamper water ice reformation. At r 30 AU, the HDO ice /H 2 O ice ratio decreases by up to one order of magnitude within 10 6 yr, without significant decrease of the water ice column density. Figure 11 shows spatial distributions of water vapor abundance and the gaseous HDO/H 2 O ratio. In our model, water vapor is moderately abundant ( 10 −8 ) in three regions: (i) the midplane inside the snowline, (ii) the disk surface at r < 50 AU, and (iii) the outer disk. They have been identified by previous disk chemical models without mixing (Woitke et al. 2009b;). This structure is conserved in the model with vertical mixing. In this subsection, we discuss how the D/H ratio of water vapor is determined in these regions.

Deuteration of Water Vapor
(i) The midplane inside the snowline (T g ∼ T d ∼ 200 K, A V > 10, and r ∼ 1 AU).
At r ∼ 1 AU, the bulk of water exists in the gas phase near the midplane, since dust temperatures are higher than the sublimation temperature of water. Since FUV photons and X-rays are strongly attenuated, H 2 O and HDO are mainly destroyed by reactions with ions and photons induced by cosmic-rays in a timescale of a few 10 5 yr. We note that cosmic-ray is also slightly attenuated; while the unattenuated cosmic ray ionization rate is 5 × 10 −17 s −1 , the ionization rate is ∼1 × 10 −17 s −1 in the midplane of r = 1 AU. H 2 O and HDO reform via proton (deuteron) transfer (e.g., H 3 O + + NH 3 ) and/or neutral-neutral reactions (OH + H 2 and OD + H 2 ; see Table 2). This cycle reduces the HDO/H 2 O ratio, since the products of HDO destruction, such as H 2 DO + , are converted not only to HDO, but also to H 2 O. The HDO/H 2 O ratio decreases from the initial value in a similar timescale to the destruction timescale, and reaches ∼10 −3 at 10 6 yr in the model without mixing. In the model with mixing, their destruction timescales are shortened by a combination of the upward transport and photoreactions at the disk surface. The HDO/H 2 O ratio decreases in a timescale of 1 × 10 4 yr (7 × 10 4 yr) in the model with α z = 10 −2 (10 −3 ), and finally reaches ∼10 −5 .
(ii) The inner disk surface (T g 300 K, A V < 1, and r 50 AU).
The HDO/H 2 O ratio is as low as or less than 10 −4 in this region, where dust and gas temperatures are decoupled. This result agrees with Willacy & Woods (2009). The HDO/H 2 O ratio is similar to the OD/OH ratio, since H 2 O and HDO is mainly formed via OH + H 2 and OD + H 2 , respectively. While the exchange reaction of OH + D deuterates OH, the reverse reaction with the activation barrier of 810 K, OD + H, prevents the significant deuteration of OH (i.e. water) at T g 300 K. Vertical mixing increases the abundance of water vapor in this region by enhancing the H 2 abundance (Section 4.1), while it does not strongly affect the HDO/H 2 O ratio.
(iii) The outer disk (T g 300 K, A V < a few, and r 50 AU).
In the region with T g < 300 K, H 2 O and HDO mainly forms via a sequence of ion-neutral reactions, followed by the recombination of H 3 O + and its isotopologues with an electron, while they are mainly destroyed via photodissociation and reactions with ions. It should be noted that the abundance of water in this region depends on the X-ray flux , and that photodesorption of water ice is important to keep a fraction of oxygen in the gas phase (e.g., Dominik et al. 2005). Vertical mixing increases both the abundances of OH and atomic deuterium by a factor of a few. The abundance of OD also increases, since it is mainly formed by the reaction of OH + D. Since H 2 DO + is formed via a sequence of ion-neutral reactions, initiated by the reactions between OD and ions, the HDO/H 2 O ratio in this region also increases in the model with mixing.
In Figure 12, we present radial variations of water vapor column density and its D/H ratio. At 3 AU r 30 AU, the column densities of water vapor increase by up to one order of magnitude via vertical mixing compared to the model without mixing. This is due to the increased abundance of hot water (T g 300 K). The column density at ∼1 AU decreases via vertical mixing by a factor of a few, since a fraction of oxygen transported from the surface is converted to O 2 . At the inner disks (< 1 AU), where gas temperature near the midplane is higher than 300 K, reformation of water vapor via Reaction (30) would be faster than O 2 formation.
Vertical mixing also affects the HDO/H 2 O column density ratio. One of the most important effects of the mixing is a substantial decrease of the ratio inside the snow line. Beyond the snow line, on the other hand, vertical mixing decreases or increases the column density ratio, depending on which dominates in the column density, hot water (T g 300 K) or warm water (T g 300 K). The ratio at r = 2 AU is high (∼10 −2 ) in the model without mixing, since a fraction of water exists in the gas phase (∼10 −6 ) near the midplane via thermal desorption of ice.
Finally, we point out that the D/H ratios of water vapor and ice are different, especially in the inner regions (r 30 AU), regardless of the strength of vertical mixing. This suggests that it is difficult to constrain the D/H ratios of water ice in the midplane of the inner disks from the observations of D/H ratios of water vapor in protoplanetary disks.

Comparisons to Other Works
The effect of mixing on water chemistry beyond the snow line (r > 10 AU) was studied by Semenov & Wiebe (2011), considering both radial and vertical mixing. Their chemical network is also based on Garrod & Herbst (2006). They found that water ice column densities in the model with and without mixing differ only by a factor of <2-5, which is consistent with our results at r 30 AU. However, in the inner disks, column densities of water ice decreases by more than one order of magnitude in our models. The difference seems to mainly come from the higher dust temperatures in our models compared to those in Semenov & Wiebe (2011). Radial mixing may also help maintain the high abundance of water ice, which will be investigated in forthcoming papers.
While there are several studies focusing on deuterium fractionation in PPDs, none of them has considered mixing. Here, we briefly compare our model without mixing to those of Willacy & Woods (2009), which focused on the inner disks (r 30 AU) considering the radial accretion. They found that both beyond and inside the snow line, the D/H ratio of water in the midplane retains its initial value (∼10 −2 ) for 10 6 yr. In our models without mixing, the D/H ratio of water ice in the midplane retains the initial value beyond the snow line, while the D/H ratio of water vapor decreases inside the snow line in a timescale of a few 10 5 yr. The different results inside the snow line seem to come from higher unattenuated cosmic-ray ionization rate in our models (5 × 10 −17 s −1 ), than that in Willacy & Woods (2009) (1.3 × 10 −17 s −1 ). The timescale of destruction and reformation of water vapor, and thus the timescale of decline of its D/H ratio is inversely proportional to the cosmic-ray ionization rate. If we adopt the unattenuated cosmic-ray ionization rate of 1.3 × 10 −17 s −1 , the timescale becomes comparable to or longer than 10 6 yr. The combination of radial accretion and evaporation of water ice at the snow line may also help maintain the initial D/H ratio of water. Willacy & Woods (2009) also found that the HDO/H 2 O ratio in the inner disk surface is as low as 10 −4 , which is consistent with our model. In order to make the molecular composition with the low HDO/H 2 O ratio, we artificially decrease the abundances of deuterated species except for HD by a factor of forty compared to the initial abundance used in our fiducial models. We reset the HD abundance so that the elemental deuterium abundance of 1.5 × 10 −5 relative to hydrogen. For species without deuterium, we use the same abundances as those in our fiducial models. Then, the elemental abundances used here are slightly different from those in our fiducial models (at most ∼4 %), except for deuterium. In this initial abundance, the HDO/H 2 O ratio is ∼5 × 10 −4 , which is similar to the cometary value (∼6 × 10 −4 ), and ∼99 % of deuterium is in HD.
In Figure 13, we show temporal variations of HDO ice /H 2 O ice ratio in the midplane of r = 50 AU. The ratio increases with time from 5 × 10 −4 to 2 × 10 −3 in 10 6 yr in the model with α z = 10 −2 , while the ratio is nearly constant in the model without mixing. Unlike the fiducial model, transport of deuterium via mixing is not important, since almost all deuterium is in HD from the beginning. Still, transport of oxygen affects water and deuterated water chemistry; atomic oxygen is transported from the surface to the deeper region, and (re)form H 2 O and HDO ices. The HDO/H 2 O ratio increases, since the D/H ratio of reformed water ice is larger than 5 × 10 −4 . After a few 10 6 yr, the HDO ice /H 2 O ice ratio reaches the steadystate value of 2 × 10 −3 , which is independent of the initial ratio. The steady-state value should correspond to the D/H ratio of reformed water ice in the disk at a given radius (i.e. temperature). Figure 14 shows radial distributions of HDO ice /H 2 O ice column density ratio at t = 10 6 yr. In the inner region (r 40 AU), where τ mix is much smaller than 10 6 yr, the HDO ice /H 2 O ice ratio is independent of the initial ratio. In the outer region, on the other hand, the ratio depends on the initial ratio at least for 10 6 yr.

Water Formation
Recently,  investigated water formation on carbonaceous grain surfaces, considering chemisorption of atomic hydrogen. They found that the reaction between physisorbed atomic oxygen and chemisorbed atomic hydrogen is very efficient at T d < 40 K, but the efficiency declines rapidly at higher temperatures. It also should be noted that chemisorption sites are unavailable on grain surfaces with ice-mantles (i.e. when x H 2 O ice > N site x d ∼ 10 −6 ). Therefore, the inclusion of water ice formation on chemisorption sites would not significantly change our results.
In our fiducial models, the desorption energy of atomic hydrogen is fixed to be 450 K, following Garrod & Herbst (2006). Watanabe et al. (2010) and Hama et al. (2012) observed various kinds of potential sites (or site distributions) with different depths (E des ) for atomic hydrogen on amorphous water ice in their experiments. According to molecular dynamics simulations, the site distribution peaks at E des ∼ 400 K and ∼600 K for crystalline and amorphous water ices, respectively (Al-Halabi & van Dishoeck 2007). Higher desorption energy of atomic hydrogen would lead to the longer residence time on grain surfaces, and thus a higher (re)formation rate of water ice. In order to check the dependence of E des on our results, we reran our model from r = 10 AU to r = 50 AU with E des = 600 K for atomic hydrogen and deuterium. The resultant radial distribution of water ice column density and its D/H ratio are shown in Figure 15. In the model with E des = 600 K, the water ice reformation compensates the destruction at r 20 AU, while in our fiducial model (E des = 450 K), it compensates the destruction only at r 40 AU. In the inner region (r 20 AU), the significant reduction of water ice column density occurs, even with E des = 600 K. The radial variations of the HDO ice /H 2 O ice column density ratio does not strongly depends on E des . If the E des is 600 K, the column density ratio of HDO ice /H 2 O ice reaches the cometary value in 10 6 yr, while the H 2 O ice abundance is kept high, in the comet forming regions (r ∼ 20-30 AU)(see Section 5.5).
In order to efficiently (re)form water ice, OH should be formed on grain surfaces as discussed in Section 4.1. In our model, OH ice is mainly formed via the reaction, The simplest reaction, O ice + H ice , is less efficient, because of the lower desorption energy (or shorter residence time on grain surfaces) of atomic hydrogen than that of HNO (Table  3). The product NO ice cycles back to HNO ice via the reaction, which is faster than O ice + H ice , because of the longer residence time of NO than that of atomic oxygen. There is a caveat about this loop of NO ice /HNO ice interconversion as pointed out by Garrod et al. (2009). The grain-surface network of Garrod & Herbst (2006), which we employed, has only a limited number of reactions involving NO ice and HNO ice . While the network would be appropriate for cold environments (10 K), where atomic hydrogen is the dominant mobile reactant, NO ice and HNO ice could be subject to various other reactions in warmer environments, where species with heavy elements also become mobile. Garrod et al. (2009) mentioned that a more comprehensive reaction network that allowed NO ice to react with, for example, atomic oxygen and nitrogen (which the current network does not) would make such loop inefficient. However, it is unclear whether the caveat applies to the disk chemistry with mixing. The physical condition of the disk, such as density, is significantly different from those of clouds/cores considered by Garrod et al. (2009). Mixing enhances the abundances of reactive species, such as atomic hydrogen and oxygen. Modeling of surface reaction on warm dust grains with a more comprehensive reaction network, which is currently unavailable, is desirable.

Photodissociation rates of ices
In our fiducial models, we calculate the photodissociation rates of ices (Equation 9), assuming that only the upper most one monolayer can be dissociated, while photoproducts immediately recombine in the deeper layers. We may underestimate the photodissociation rates of ices, if the probability of recombination of photofragments in the deeper layers, which would depend on ice temperatures (Öberg et al. 2009c), is much less than unity. To investigate the dependence of our results on the parameter N p , we rerun the calculation in which N p is set to be equal to the number of monolayers (N layer ).
The increased photodissociation rate of water ice enhances the conversion of water ice to CO 2 ice, since OH ice produced by water ice photodissociation are partly converted to CO 2 ice. Then, the region where CO 2 ice becomes the dominant oxygen reservoir is extended to r ∼ 50 AU in the model with α z = 10 −2 , while the boundary was r ∼ 40 AU in the fiducial model. In our model, CO 2 ice is mainly formed by the reaction of CO ice + OH ice . The size of the CO 2 ice dominated region would also depend on the value of the activation energy barrier of this reaction. The barrier is assumed to be 80 K (Ruffle & Herbst 2001) in our models, while Noble et al. (2011) concluded that the reaction is likely to have a higher barrier based on their experiments. If we adopt a higher value, the CO 2 ice dominated region should become smaller. In our fiducial model with N p = 2, we artificially switch off the D ice + OH ice branch of HDO ice photodissociation (see Section 3.2.1). When we assume N p = N layer , there is no reason to switch off the branch; here, we assume that the branches to produce OH ice and OD ice are equally weighted. We found that the HDO ice /H 2 O ice ratios at t = 10 6 yr in the model with α z = 10 −2 are smaller than those in our fiducial model by up to a factor of five in the regions where water ice is the dominant oxygen reservoir (r 50 AU).

Effect of Accretion and Grain Evolution
In the present study, we consider only vertical mixing as the mass transport in PPDs. However, radial mass transport (in the combination with the vertical mixing) could be also important for disk chemistry (e.g., Tscharnuter & Gail 2007;Nomura et al. 2009;Heinzeller et al. 2011;Semenov & Wiebe 2011). In Figure 16, we compare τ mix with the timescale of radial accretion in the model with α z = 10 −2 . The accretion timescale is estimated by where Σ is the surface density of gases. Note that when we determine the radial distribution of surface density in our disk models, we assume a viscous parameter of α = 10 −2 as mentioned in Section 2. Figure 16 shows that the accretion timescale is comparable to τ mix . The radial accretion would increase the water abundance and the HDO/H 2 O ratio in our fiducial model with mixing at r 30 AU, since they are higher at the outer radii.
We also neglect grain evolution in this work. In our chemical models, it is assumed that 0.1 µm dust grains are uniformly distributed in the disk with the dust-to-gas mass ratio of 10 −2 . In reality, PPDs contain grains as large as 1 mm, and a fraction of them would be settled to the midplane (e.g., Przygodda et al. 2003;Furlan et al. 2006Furlan et al. , 2011. While the 0.1 µm grains are well-coupled to the gas, the coupling is less efficient for larger grains. If large grains with ice mantles settle the midplane, but small grains remain in the disk surface to shield UV radiation, destruction of water ice via the combination of vertical mixing and photoreactions would become less efficient, although water vapor is subject to the destruction as long as turbulence exists. A similar situation is considered to account for the weak emission line of water vapor detected towards TW Hya (Hartogh et al. 2011).

Effect of Dead Zone
In our fiducial models, we assume that α z is constant in space. According to MHD simulations, however, velocity dispersion in a dead zone is smaller than that in a MRI active region by more than one order of magnitude, depending on the strength of the magnetic field (Okuzumi & Hirose 2011;Gressel et al. 2012). In this subsection, we consider the models, in which α z is dropped by one order of magnitude in a dead zone. We define a dead zone as the regions with the magnetic Reynolds number, R M , is less than 100 (e.g., Fleming et al. 2000). Following Perez-Becker & Chiang (2011), the magnetic Reynolds number is evaluated by where h, η, and x ′ e are the pressure scale height, electrical resistivity of the gas, and the electron abundance with respect to neutrals, respectively. For example, the threshold values of x ′ e for the dead zone in our disk model are 9 × 10 −13 and 4 × 10 −13 at r = 5 AU and 10 AU, respectively.
In Figure 17, we show a dead zone boundary (solid line), plotted over the distribution of water ice abundance in the model with α z = 10 −2 at t = 10 4 yr. The dead zone exists near the midplane (z/r 0.1) at r 10 AU. The dashed line in Figure 17 indicates the height z * , at which the effective flux of water ice to the disk surface is determined (Section 4.2). The height z * is located above the dead zone boundary. In other words, the water ice layer is thicker than the dead zone, and the net upward transport of water ice takes place in the MRI active region in our models. We performed a calculation considering the dead zone, and confirmed that the existence of the dead zone does not significantly change our results on both water ice and vapor.

Cometary Water
Comets observed today are believed to be supplied from two distinct reservoirs, the Oort cloud and the Kuiper belt. Until recently, it was widely accepted that the Oort cloud comets (OCCs) formed in and was scattered outward from the giant planet forming region (r ∼ 5-30 AU), while Jupiter-family comets (JFCs) formed in the Kuiper belt region (r 30 AU) (Brownlee 2003). However, recent dynamical models of solar system evolution suggest more complicated scenarios; substantial migration of the giant planets occurred and it led to large scale mixing of distributions of planetesimals (Gomes et al. 2005;Walsh et al. 2011). If this is the case, OCCs and JFCs would share their origins at least in part.
In our fiducial models with α z = 0 and 10 −3 , in which the initial HDO/H 2 O ratio is ∼2×10 −2 , the model D/H ratios in water ice remain much higher than the cometary value even at t = 10 6 yr. In the model with α z = 10 −2 and E des = 600 K for atomic hydrogen, water ice is abundant (∼10 −4 ) and its D/H ratio is comparable to the cometary value at r ∼ 20-30 AU as shown in Figure 15. Such region with abundant water ice and the preferable D/H ratio can extend to inner radius, if the reformation of water ice is more efficient than that in our models. Then, our model suggests the possibility that the D/H ratio of cometary water could be established (i.e. cometary water could be formed) in the solar nebula, even if the D/H ratio of water ice formed in the parent molecular cloud/core was very high (∼10 −2 ) as observed in NGC 1333-IRAS2A and NGC 1333-IRAS4A (Taquet et al. 2013a This contradicts with the classical dynamical model in which the OCCs originated in the inner region than the JFCs, and may support the recent dynamical model with migration of giant planets. It should be noted, however, that the radial gradient of HDO/H 2 O ratio could be reversed, if grain settling proceeds faster in inner radii. In addition, the radial distribution of the HDO/H 2 O ratio can also be affected by the radial transport of low D/H water near the central star, which is not yet considered in our model (Yang et al. 2013). Further studies are needed to constrain the formation regions of JHCs and OCCs from their HDO/H 2 O ratios.

CONCLUSION
We have investigated water and deuterated water chemistry in protoplanetary disks irradiated by UV and X-ray from a central T Tauri star. We have solved chemical rate equations with the diffusion term, mimicking the turbulent mixing in the vertical direction. Oxygen is mainly in atomic form in the disk atmosphere, while it is in water near the midplane. When turbulence exists, water near the midplane is transported to the disk surface and destroyed by photoreactions, while atomic oxygen is transported to the midplane and reforms water and/or other molecules. We found that this cycle significantly affects water and deuterated water chemistry. Our conclusions are as follows.
1. Beyond the snow line, the cycle decreases the column densities of water ice by more than one order of magnitude within 10 6 yr in the inner disk (r 30 AU), where dust temperatures are too high to form OH radical on grain surfaces. Once OH radical is formed in the gas phase, it is converted to O 2 via the reaction of OH + O before it is adsorbed onto dust grains to form water ice, because of the high abundance of atomic oxygen near the O/H 2 O ice transition. The outer edge of such regions moves to r = 20 AU, if the desorption energy of atomic hydrogen is as high as 600 K. Our model indicates that water ice could be deficient even outside the sublimation radius.
2. At r 30 AU, the cycle decreases the D/H ratios in water ice from ∼2×10 −2 , which set by the collapsing core model, to 10 −3 -10 −2 within 10 6 yr, without significant decrease of the column densities. The resultant ratio depends on the strength of mixing and the radial distance from the central star. If E des is 600 K for atomic hydrogen, the D/H ratio of water ice decreases to 10 −4 -10 −3 at r ∼ 20-30 AU without significantly decreasing the water ice column density. Our model suggests that the D/H ratio of cometary water could be established (i.e. cometary water could be formed) in the solar nebula, even if the D/H ratio of water ice formed in the parent molecular cloud/core was very high (∼10 −2 ) as observed in some hot corinos.
3. We confirmed that water vapor has moderately high abundance ( 10 −8 ) in the three regions: (i) the midplane inside the snowline, (ii) the inner disk surface, and (iii) the outer disk as shown in previous works. We found that this structure conserves in the model with vertical mixing. Inside the snow line, the D/H ratios in water vapor become as low as 10 −5 in the model with mixing.     -Spatial distributions of water ice abundance in the model with α z = 0 (left panels), 10 −3 (middle panels), and 10 −2 (right panels) at 10 5 yr (top panels) and 10 6 yr (bottom panels). The vertical axes represent height normalized by the radius. The dashed lines indicate vertical visual extinction of unity and ten.       11.-Spatial distributions of water vapor abundance (blue colors) and its D/H ratio (orange colors) in the model with α z = 0 (left panels), 10 −3 (middle panels), and 10 −2 (right panels) at 10 6 yr. The dashed lines indicate vertical visual extinction of unity and ten. The D/H ratios are only shown in regions where the water vapor abundance is higher than 10 −12 .