11institutetext: Faculty of Aerospace Engineering, Delft University of
Technology, Delft, The Netherlands 22institutetext: School of Physics &
Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK
33institutetext: Kapteyn Astronomical Institute, University of Groningen, P.O.
Box 800, 9700 AV Groningen, The Netherlands;
33email<EMAIL_ADDRESS>44institutetext: Max Planck Institute for
Extraterrestrial Physics, Giessenbachstrasse 1, 85741 Garching, Germany
55institutetext: Leiden Observatory, Leiden University, P.O. Box 9513, NL 2300
RA Leiden, The Netherlands 66institutetext: Centre for Exoplanet Science,
University of St Andrews, North Haugh, St Andrews, KY16 9SS, UK
77institutetext: Universitäts-Sternwarte, Fakultät für Physik, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany

# Ices in planet-forming disks: Self-consistent ice opacities in disk models

Aditya M. Arabhavi 11223366 Peter Woitke 2266 Stéphanie M. Cazaux 1155 Inga
Kamp 33 Christian Rab 4477 Wing-Fai Thi 44

(Received July 21, 2021; accepted August 18, 2022)

###### Abstract

Context. In cold and shielded environments, molecules freeze out on dust grain
surfaces to form ices such as $\rm H_{2}O$, $\rm CO$, $\rm CO_{2}$, $\rm
CH_{4}$, $\rm CH_{3}OH$, and $\rm NH_{3}$. In protoplanetary disks, such
conditions are present in the midplane regions beyond the snowline, but the
exact radial and vertical extension depend on disk mass, geometry, and stellar
ultra-violet irradiation.

Aims. The goal of this work is to present a computationally efficient method
to compute ice and bare-grain opacities in protoplanetary disk models
consistently with the chemistry and to investigate the effect of ice opacities
on the physico-chemical state and optical appearance of the disk.

Methods. A matrix of Mie efficiencies is pre-calculated for different ice
species and thicknesses, from which the position dependent opacities of icy
grains are then interpolated. This is implemented in the PRODIMO code by a
self-consistent solution of ice opacities and the local composition of ices,
which are obtained from our chemical network.

Results. Locally, the opacity can change significantly, for example, an
increase by a factor of more than 200 in the midplane, especially at ultra-
violet and optical wavelengths, due to ice formation. This is mainly due to
changes in the size distribution of dust grains resulting from ice formation.
However, since the opacity only changes in the optically thick regions of the
disk, the thermal disk structure does not change significantly. For the same
reason, the spectral energy distributions (SEDs) computed with our disk models
with ice opacities generally show only faint ice emission features at far-IR
wavelengths. The ice absorption features are only seen in the edge-on
orientation. The assumption made on how the ice is distributed across the
grain size distribution (ice power law) influences the far-infrared and
millimeter slope of the SED. The ice features and their strengths are
influenced by the ice power law and the type of chemistry. Our models predict
stronger ice features for observations that can spatially resolve the disk,
particularly in absorption.

###### Key Words.:

protoplanetary disks – opacity – methods: numerical

##  1 Introduction

Protoplanetary disks are made up of gas and dust with different spatial
distribution and temperature structures. The midplane is generally cold and
dense for both dust and gas components (e.g. Bergin et al., 2007; Öberg &
Bergin, 2021). In those cold regions, which are well shielded from both
stellar and interstellar ultra-violet (UV) radiation fields, gas phase species
can accrete onto dust grain surfaces to form ices. The local thermo-chemical
conditions in these regions determine the volume and chemical composition of
the forming ice. This ice accretion affects the grain size and composition
which are important for opacity calculations (Ossenkopf & Henning, 1994;
Allamandola et al., 1999; Boogert & Ehrenfreund, 2004; Öberg, 2009).

Observational data of dust in the diffuse interstellar medium, molecular
clouds, and young stellar objects show the presence of ice covering the dust
grains (Merrill et al., 1976; Whittet et al., 1983, 1998; Chiar et al., 2000,
2002; Boogert et al., 2008; Goto et al., 2018). For protoplanetary disks, the
resolution, spectral range, and sensitivity requirements limit the
observational capabilities. Yet some observations of $\rm H_{2}O$ ice in
absorption at 3$\mu$m (for example Terada et al., 2007; Honda et al., 2009,
2016; Terada & Tokunaga, 2017) and emission at 47 $\mu$m and 63 $\mu$m
(McClure et al., 2012, 2015; Min et al., 2016a) are available. Thus, disk
modeling becomes important for understanding the properties of icy dust in
disks as well as to design observing strategies for observing ices.

Observations, such as Terada et al. (2007) and Terada & Tokunaga (2012),
indicate that the observed profile of an ice feature (for example the 3 $\mu$m
$\rm H_{2}O$ ice feature) is a combination of diverse factors. In addition to
the contributions of the ice species themselves, which include presence of
different ice species such as $\rm CO$, $\rm CH_{4}$, $\rm NH_{3}$, $\rm
CH_{3}OH$ (Brooke et al., 1999; Dartois & d’Hendecourt, 2001; Dartois et al.,
2002), optical behavior of large grains (Smith et al., 1989; Boogert et al.,
2015) and thermal state of the ices (Boogert et al., 2015; Sandford &
Allamandola, 1988; Hudgins et al., 1993; Palumbo & Strazzulla, 1993)
contribute to the shape of the ice feature. Grain growth is also expected in
cold and dense environments (e.g. Terada & Tokunaga, 2012; Terada et al.,
2012; McClure et al., 2012).

Generally, disk models (e.g. Nomura & Millar, 2005; D’Alessio et al., 2006;
Semenov & Wiebe, 2011; Bruderer et al., 2012; Woitke et al., 2016) assume a
grain size distribution with a minimum and maximum grain size, and a power law
with index $a_{\rm pow}$ that is either obtained by fitting the observations
or set to 3.5 (MRN-distribution, Mathis et al., 1977). These grains are
considered to be composed of refractory materials such as silicates and
conducting materials. Some models consider the presence of ices by using ice
opacities (e.g. Nomura & Millar, 2005; Pontoppidan et al., 2005; McClure et
al., 2012; Kamp et al., 2018). However, some of these models (Nomura & Millar,
2005; D’Alessio et al., 2006; McClure et al., 2012, 2015) consider segregated
grain distribution, that is, having separate size distribution of ice grains
apart from refractory grains. Moreover, for calculating the opacities, these
models (Nomura & Millar, 2005; D’Alessio et al., 2006; McClure et al., 2015)
either assume that ice exists everywhere below sublimation temperature,
neglecting photodesorption, or do not consistently solve for photodesorption
and ice abundance. Further, Nomura & Millar (2005), D’Alessio et al. (2006),
McClure et al. (2012), McClure et al. (2015), and Kamp et al. (2018) consider
only $\rm H_{2}O$ ice and Pontoppidan et al. (2005) consider $\rm H_{2}O$,
$\rm CO$, and $\rm CO_{2}$ ice species for opacity calculations. Segregated
grain distribution do not account for increase in grain size due to ice
accretion. Considering only $\rm H_{2}O$ ice underestimates the effect of ice
accretion on size distribution of grains. A more recent study of ice opacities
in disk models is presented in Ballering et al. (2021), where six ice species
($\rm H_{2}O$, $\rm CO$, $\rm CO_{2}$, $\rm CH_{3}OH$, $\rm NH_{3}$, and $\rm
CH_{4}$) are considered along with two populations of refractory grains. The
disk region is spatially classified into 25 zones which can be approximated
with the same opacities. Thus avoiding the computationally intensive, yet more
accurate, point-by-point opacity calculation. The effects of an interstellar
UV field is also neglected. The model produces substantial ice absorption
features including the 3 $\mu$m water ice band in disks that are not edge on.
A more consistent treatment of the effect of ice on dust composition and size,
hence, will help improve the understanding of ice observations in disks.

Ossenkopf & Henning (1994) show, for a fixed volume of water ice, that the
opacity changes when ice mantles are present on dust grains, which is not
taken into consideration by many disk models. Kamp et al. (2018) assumes that
$\rm H_{2}O$ ice is present in region defined by the snowline (Min et al.,
2016a) which takes into account the temperature, pressure, and UV radiation
field. The ice is distributed as ice mantles on top of underlying bare grain
size distribution. Ices are assumed to be a constant mass fraction of the
total dust mass, thus the ice abundance is not informed by the local chemical
state. This leads to having too much ice in the upper midplane or too little
ice in the inner midplane. Pontoppidan et al. (2005) use a single refractory
grain size population, upon which ice mantle composed of $\rm H_{2}O$, $\rm
CO$, and $\rm CO_{2}$ ices is present. However, the model considers that ices
evaporate at 90K (20K for CO) and also the local ice abundances are obtained
by ice mixtures of fixed proportions ($\rm H_{2}O$:$\rm CO_{2}$:$\rm CO$=$\rm
1$:$\rm 0.2$:$\rm 0.03$ and inclusions of $\rm CO_{2}$:$\rm CO$=$\rm 1$:$\rm
0.7$). In this paper, we present a model in which, for the first time, the ice
chemistry, opacities, and radiative transfer are treated self-consistently.

We present here a computationally efficient method of including position
dependent opacities in disk models, to account for ice abundances based on
local thermo-chemical state, and discuss the implications for observations. In
Sect. 2 the description of the disk model, PRODIMO, is provided in which the
position dependent ice opacities are implemented. In Sect. 3, the method of
calculating position dependent opacities is explained along with the
verification of the calculation. The effect of the position dependent opacity
on a typical T Tauri disk is discussed in Sect. 4. We also address two main
factors that affect the ice features in our disk models, namely the relation
between ice thickness and grain size, and the chemistry. We discuss the
results of our models in the context of ice observations in Sect. 5. The
method and results are summarized and concluded in Sect. 6.

##  2 Disk model and optical data

PRODIMO111Version v2.0.0, https://prodimo.iwf.oeaw.ac.at/ or Protoplanetary
Disk Model (Woitke et al., 2009; Kamp et al., 2010; Thi et al., 2011; Rab et
al., 2018) is used in this paper to implement the position-dependent dust and
ice opacities. The existing model, as explained in Woitke et al. (2009) and
Woitke et al. (2016), uses a fixed grain size distribution with no ices to
calculate the opacities. The model starts by setting up the physical state of
the disk such as gas and dust densities including dust settling of the bare
grains. The code then calculates the dust opacities and solves the continuum
radiative transfer, thereby determining the dust temperature structure.
Subsequently, the code computes the gas phase and ice chemistry and solves the
heating & cooling balance to find the gas temperature structure. We use steady
state chemistry in all the models presented in this paper. One model also
includes grain surface reactions. More details about the chemistry and the
disk parameters are given in Appendix B.

The bare grain opacities are explained in Woitke et al. (2016), based on Min
et al. (2016b). The opacity calculations assume that the grains are spherical
and the materials are well mixed. The refractory dust grains are assumed to be
composed of 60% amorphous silicate (Dorschner et al. (1995), $\rm
Mg_{0.7}Fe_{0.3}SiO_{3}$) and 15% amorphous carbon (Zubko et al., 1996), with
a porosity of 25%.

The computation method presented in this paper requires multiple iterations to
converge. The first iteration uses bare grain opacities. The resulting
position-dependent ice abundance and position-dependent settled bare grain
size distribution are used to recompute the dust opacities including the ices
in the next iteration. Further, the code continues to iterate between
radiative transfer, chemistry with ice formation, and opacity calculation,
until convergence is achieved. The models presented in this paper converged
within 3 such iterations.

The opacity calculations are based on Mie theory as explained in Woitke et al.
(2016). The dust grains (including those coated with ice) are considered to be
spherical. Further, Woitke et al. (2016) use Mie theory and distributed hollow
spheres (Min et al., 2005) for calculating opacities. However, the grain size
distribution at each grid point in the disk changes and the distributed hollow
spheres are computationally expensive, in this paper distributed hollow
spheres are not used. In this paper, the dependence of ice thickness on grain
size is generalized and made computationally easy to apply to disk models. The
porosity of the refractory grains and ice layer is fixed to a certain constant
value, that is, the possible enhancement or diminution of porosity upon ice
accretion, grain collision, shattering, and other thermal events is neglected.

##  3 Opacity: Ice composition, thickness, and refractory materials

Extinction opacity ($\kappa_{\lambda}^{\rm ext}$) can be defined as the total
absorption ($\kappa_{\lambda}^{\rm abs}$) and scattering
($\kappa_{\lambda}^{\rm sca}$) cross sections of the dust distribution, as
shown in Eq. (1) (Krugel, 2002):

$\displaystyle\kappa_{\lambda}^{\rm ext}$ $\displaystyle=\kappa_{\lambda}^{\rm
abs}+\kappa_{\lambda}^{\rm sca}$ (1) $\displaystyle Q_{\lambda}^{\rm ext}$
$\displaystyle=Q_{\lambda}^{\rm abs}+Q_{\lambda}^{\rm sca}$
$\displaystyle\kappa_{\lambda}^{\rm ext}$ $\displaystyle=\int_{a_{\rm
min}}^{a_{\rm max}}\pi a^{2}Q^{\rm ext}(a,\lambda)f(a)da$

where, $Q^{\rm ext}(a,\lambda)$, $Q^{\rm abs}(a,\lambda)$, $Q^{\rm
sca}(a,\lambda)$ are the extinction, absorption, and scattering efficiencies
defined as the ratio of respective cross sections to the geometric cross
section of the grain, $a$ is the grain size, $a_{\rm min}$ and $a_{\rm max}$
are the minimum and maximum grain sizes, $f(a)$ is the considered size
distribution ($=\frac{dN}{da}$, i.e., number of grains, $dN$, per size bin
$da$).

Due to ice accumulating on the bare (refractory) grains, both extinction
efficiency $Q(a,\lambda)$ and grain size distribution $f(a)$ change. These
directly affect the opacity $\kappa_{\lambda}$. The effect of ice on these two
factors are discussed in the following subsections.

###  3.1 Effect of ice on size distribution

The abundance of ice varies throughout the disk (see appendix). This means
that the thickness of ice layer differs at each point in the disk. The
thickness of ice on the grains at any point depends on the local volume of ice
available and the total grain surface area available for ice accretion at that
point in the disk. A simple relation can be written as shown in Eq. (2).

$\displaystyle V_{\rm IG}$ $\displaystyle=V_{\rm I}+V_{\rm BG}$ (2)
$\displaystyle n_{\rm d}\int_{a_{\rm min}}^{a_{\rm
max}}\frac{4}{3}\pi(a+\Delta a)^{3}f(a)da$ $\displaystyle=V_{\rm I}+n_{\rm
d}\int_{a_{\rm min}}^{a_{\rm max}}\frac{4}{3}\pi a^{3}f(a)da$

where, $V$ is volume of material per unit volume of space with subscripts $\rm
IG$ for icy grains, $\rm I$ for ice, and $\rm BG$ for bare grains, $a$ is the
bare grain size, $\Delta a$ is the ice layer thickness, $f(a)$ [$\rm cm^{-1}$]
is the normalized size distribution function of the grains, $a_{\rm min}$ and
$a_{\rm max}$ are the minimum and maximum sizes of the bare grains, and
$n_{\rm d}$ [$\rm cm^{-3}$] is the number density of grains. We consider the
following two cases. (1) If, in Eq. (2), $\Delta a$ is constant for all grain
sizes, then the ice/refractory volume ratio ($V_{\rm I}$/$V_{\rm BG}$) is
significantly larger for small grains than for large grains. (2) Ice thickness
can also depend on the grain size. This can be due to coagulation or
shattering of icy grains in the midplane or heating and re-condensation events
(the latter leads to growth of ice faster on larger grains, Kuroiwa & Sirono,
2011). In this paper we assume that the limiting case of such scenario is an
uniform ice/refractory volume ratio for grains of different sizes.

Figure 1: Icy grain size distribution for limiting values of ice coating
powerlaw, $q$. Bare size distribution: $a_{\rm min}$ = 0.05 $\mu$m, $a_{\rm
max}$ = 3mm, powerlaw index 3.5, $V_{\rm BG}=6.2894\times
10^{-13}\rm\,cm^{3}/cm^{3}$, $V_{\rm I}=2.7577\times
10^{-12}\rm\,cm^{3}/cm^{3}$.

In order to allow for a general approach, which includes both limiting cases,
we assume that the resulting grains are still spherical. For an icy grain of
certain refractory size $a$, if $V_{\rm I}$, $V_{\rm BG}$, and $V_{\rm IG}$
are volumes of ice, refractory material, and total volume of the icy grain,
respectively, then the ratio ${V_{\rm I}}/{V_{\rm IG}}$ and ${V_{\rm
BG}}/{V_{\rm IG}}$ are the volume fractions of ice and refractory material in
an icy grain respectively. If a grain with refractory size $a$ is present in
an environment where there is no coagulation, shattering, or heating and
recondensation events, then the ice thickness $\Delta a$ follows some constant
$t$ and hence the volume fractions change across the size distribution.
However, if the grain is present in an environment where the grain size
distribution has been reached by coagulation, shattering, or heating and
recondensation events, then the limiting case (according to the assumption)
leads to an ice thickness $\Delta a$ proportional to the size of refractory
material, $t\cdot a$, thus preserving the ice/refractory volume fraction. The
following powerlaw approach is used to include both limiting cases

$\Delta a=t\cdot a^{q}$ (3)

where, $t$ and $q$ are arbitrary coefficients that are indicative of local ice
abundance and the local dominance of shattering, coagulation, heating and
recondensation events respectively. The two limiting cases are (1)
$q\\!=\\!0$, where the thickness of the ice layer is constant and (2)
$q\\!=\\!1$, where the ice/refractory volume fraction is constant and hence
the thickness linearly depends on the bare grain size. The values $0<q<1$
relate to cases where the grain size distribution is affected by coagulation,
shattering, or heating and recondensation events in varying levels. We assume
that $f(a)$ in Eq. (2) is the size distribution of the underlying bare grains.

The effect of the ice covering the dust on the size distribution (dust + ice)
is shown in Fig. 1. For a constant ice thickness ($q\\!=\\!0$), the relative
increase in size of the smaller grains is several orders of magnitude and the
new minimum grain size is determined by the ice layer thickness. The number
density of grains of size close to the ice thickness increases significantly.
However, the change in size for larger grain sizes is negligible. The overall
effect on the grain size distribution is that it does no longer follows a
powerlaw. For the other limiting case $q\\!=\\!1$, the small grains have small
ice thicknesses and large grains have large ice thicknesses. Both minimum and
maximum grain sizes are changed, but this is not as significant as in the case
of constant ice thickness. Further, the total grain surface area increases,
still dominated by smaller grains, but the size distribution powerlaw is
retained.

The total volume of ice at any point in the disk $V_{I}$ is obtained from the
chemical output of the disk code. If $n_{\rm i}\rm\,[cm^{-3}]$ is the number
density at a grid point, $\mu_{\rm i}\rm\,[g]$ is the molecular mass,
$\rho_{i}\rm\,[g/cm^{3}]$ is the material mass density of ice species $\rm i$,
and $V_{\rm i}$ is the volume of ice species $\rm i$, then the total volume of
ice per unit volume of space, $V_{\rm I}$, at a grid point is given by Eq. (4)
for N ice species.

$\displaystyle V_{\rm i}=\frac{n_{\rm i}\mu_{\rm i}}{\rho_{\rm i}}$ (4)
$\displaystyle V_{\rm I}=\sum_{{\rm i}=1}^{N}V_{\rm i}$

Using Eq. (4), Eq. (2) is solved numerically for $t$ at each grid point in the
disk using Newton-Raphson method by approximating the root of the equation
$V_{\rm IG}-V_{\rm I}-V_{\rm BG}=0$, for a given ice power law, $q$.

###  3.2 Effect of ice on extinction, absorption, and scattering efficiency

Figure 2: Extinction, absorption, and scattering opacity of bare grains
(dashed) and ice coated grains (solid). A constant thickness of ice layers is
assumed for all grains ($q\\!=\\!0$). The composition of grains is $\rm
Mg_{0.7}Fe_{0.3}SiO_{3}$, amorphous carbon, and porosity in the ratio
0.60:0.15:0.25, and that of ices: $\rm H_{2}O$, $\rm CH_{4}$, $\rm NH_{3}$ in
the ratio 0.60:0.20:0.20. The volume of dust is $V_{\rm BG}=6.2894\times
10^{-13}\rm\,cm^{3}/cm^{3}$ and thickness of ice is $\Delta a=0.246\,\mu m$.
Refractory dust grain size distribution parameters are $a_{\rm
min}\\!=\\!0.05\,\mu$m, $a_{\rm max}\\!=\\!3$ mm, and powerlaw index 3.5.
These opacities have been calculated using effective medium theory.

Extinction ($Q_{\lambda}^{\rm ext}$), scattering ($Q_{\lambda}^{\rm sca}$),
and absorption ($Q_{\lambda}^{\rm abs}$) efficiencies are related as given in
Eq. (1). These efficiencies are complex functions of size parameter ($x=2\pi
a/\lambda$) and complex refractive index $m$($\lambda$).

The default Mie-efficiency calculations in PRODIMO, as explained in Sect. 2,
assume that the same grain material is present everywhere in the disk (i.e.
position-independent). In this work, however, at all locations in the disk
where the local physical state allows for the presence of ice on grains, the
opacity calculations include the increase of size and change of material
composition of the grains due to ice formation. Two factors that should be
considered for efficiency calculations are the different ice species and their
abundances. The former affects the complex refractive index and the latter
affects both the size and the refractive index of the icy grains. Figure 2
shows the effect of these two factors on extinction, absorption, and
scattering opacities by assuming a grain size distribution and fixed ice
thickness ($q\\!=\\!0$) and composition. The figure shows the opacity cross
sections per gram of refractory material [cm2/g]. Scattering opacity increases
both at shorter and longer wavelengths and absorption increases in the mid- to
far-infrared range. This leads to an extreme albedo of ¿99% at the optical to
near-infrared wavelength region. Estimating the effects of these two factors,
the different ice species and their abundances, on the Mie efficiencies at
each grid point are individually explained in the following subsections.

####  3.2.1 Effect of ice thickness on Mie efficiencies

Figure 3: Extinction efficiency for a grain size 0.0559$\mu$m (with respect to
the bare grain geometrical cross section) with different ice thicknesses
composed of $\rm H_{2}O$, $\rm NH_{3}$, $\rm CO$, $\rm CO_{2}$, $\rm CH_{4}$,
$\rm CH_{3}OH$ in ratio 0.45:0.20:0.10:0.05:0.05:0.15 with 20% of porosity.
Also shown is the validity of the interpolation for a constant ice thickness
layer of 0.515 $\mu$m. The embedded panel shows the relative error between the
interpolation and explicit calculation of efficiency for the icy grain.

To estimate the position dependent opacities, a Mie efficiency matrix
($Q$-matrix) is first generated. Each point in this multi-dimensional matrix
corresponds toMie efficiency defined by individual ice species, ice thickness,
ice power law, bare grain size, and wavelength. For a given bare grain size
distribution and the wavlength bins, a finite number of ice thickness and ice
power law values (Eq. 3) for each ice species are used to populate this
matrix. Figure 3 shows how extinction efficiency changes as the thickness of
ice layer increases. Since it is computationally expensive to calculate the
efficiencies based on thickness at each grid point using Mie routines, log
interpolation (Eq. 5) between points in the pre-calculated $Q$-matrix is done.

In Eq. (5), $\theta$ and $\phi$ are the interpolation weights which define the
closeness of the local thickness to the nearest two thickness grid points.
However, for interpolation between the grid points of the smallest ice
thickness and bare grain opacities, a simple linear interpolation is used.
Figure 3 presents the values obtained by such interpolation. It should be
noted that in Fig. 3, the efficiency, $Q_{\rm ext}$, is the ratio of
extinction cross section of ice covered grains to the geometrical cross
section of bare grains. Here, bare grain geometrical cross section is
considered because it makes the log interpolation of efficiency simple.

$\displaystyle Q_{\rm a,\lambda}(r,z)=Q_{\rm a,\lambda}(t_{1})^{\theta}\cdot
Q_{\rm a,\lambda}(t_{2})^{\phi}$ (5) $\displaystyle\theta=\frac{log(t_{\rm
2}/t_{\rm r,z})}{log(t_{\rm 1}/t_{\rm 2})}\ \ \ \phi=\frac{log(t_{\rm
1}/t_{\rm r,z})}{log(t_{\rm 1}/t_{\rm 2})}$

where, $\rm r,z$ are the radial distance from the rotation axis and vertical
distance from the midplane of a point in disk, $t_{\rm r,z}$ is the thickness
coefficient (see Eq. 3) of ice layer at the location $\rm(r,z)$ in the disk
for which the opacity is being estimated, $t_{\rm 1}$ and $t_{\rm 2}$ are the
two ice thickness values considered in the pre-calculated $Q$-matrix closest
to $t_{\rm r,z}$, using which the efficiencies $Q_{\rm a,\lambda}(t_{\rm 1})$
and $Q_{\rm a,\lambda}(t_{\rm 2})$ are obtained from the matrix, and $\theta$,
$\phi$ are the weights for log interpolation. Finally, the values calculated
from Eq. (5) are scaled by $(a$/$\bar{a})^{2}$ to obtain the actual
efficiencies with respect to ice covered grains ($\bar{a}$ is size of the
grain with ice. i.e. $\bar{a}=a+\Delta a=a+t_{\rm r,z}\cdot a^{q_{\rm r,z}}$).
However, the local ice thickness power law ($q_{\rm r,z}$) requires extensive
solutions of shattering, coagulation physics, and heating-recondensation
events, which cannot be obtained by disk chemistry. Hence, trivial values
($q$=0,1) are considered in Section 4 independent of position (i.e. same
throughout the disk).

####  3.2.2 Effect of ice composition on Mie efficiencies

The protoplanetary disks are generally dominated by $\rm H_{2}O$ ice (Öberg et
al., 2011). Though the direct observation of ices in disks are scarce,
observations of molecular clouds and young stellar objects (YSOs) show the
presence of $\rm H_{2}O$, $\rm CO$, $\rm CO_{2}$, $\rm NH_{3}$, $\rm
CH_{3}OH$, $\rm CH_{4}$ amongst other ice species (Gibb et al., 2000, 2004;
Boogert et al., 2008; Pontoppidan et al., 2008; Noble et al., 2013). These six
ice species are considered for opacity calculations in this paper. Our sources
for the refractive indices are listed in Table 1. In this paper, the
composition of ice at any point in the disk (the number density of ice
species, see Eq. 2) is obtained from the chemistry results of the disk code,
which is explained in Kamp et al. (2017). The ice composition is determined by
accretion and desorption processes (including thermal, UV, and cosmic ray
processes). The effect of grain surface chemistry is discussed in Section 4.

It is found that the efficiencies based on local ice composition, for a given
ice thickness and a given bare grain size distribution, can be approximated by
a simple volume fraction weighted summation of efficiencies from the
$Q$-matrix for grains coated with different pure ices. This is shown in Eq.
(6). This method is computationally efficient because the number of Mie
calculations required becomes independent of the spatial grid size of the
disk, instead depends on the number of ice species, $N$, considered.

$Q_{\rm a,\lambda}({\rm
Bare,Ice_{1},Ice_{2},...,Ice}_{n})=\sum_{i=1}^{N}\left(\eta_{i}Q_{\rm
a,\lambda}({\rm Bare+Ice}_{i})\right)$ (6)

where, $Q_{\rm a,\lambda}$ is the efficiency, the left hand side of the
equation is the efficiency of the grains coated with local mixed ice, composed
of species $1$ to $N$, $\eta_{i}$ is their respective ice abundance volume
fraction ($\eta_{\rm i}=V_{\rm i}/V_{\rm I}$, Eq. 4). It should be noted that
in Eq. (6), the $Q_{a,\lambda}$ on left and right side of the equation are
based on the same refractory grain composition and size. Further,
$Q_{a,\lambda}$ on the right is calculated using the same thickness of ice as
in the left except that the ice layer is assumed to be entirely composed of
$\rm Ice_{i}$. The validity of this approach is shown in Fig. 4, which
compares the efficiencies calculated for grain coated with a mixed ice layer
(1) by using effective medium theory for ice mixture and (2) by weighted
summation as given in Eq. (6). The relative error plotted in the bottom panel
of the figure shows that the interpolation is in good agreement with
calculation using effective medium Mie routine. These approximations (Eqs. 5
and 6) introduces uncertainties (of similar order as the relative errors) in
the ice features in the disk spectra obtained from the solution. However,
these uncertainties are $\sim$1 order of magnitude smaller than the
uncertainties introduced due to difference in the opacities of a mixed and a
differentiated core-mantle grain. Hence, these approximations are sufficiently
accurate for our models but the uncertainties have to be taken into account
while fitting observations. The following section discusses the effective
medium and core-mantle Mie theory.

Figure 4: Extinction efficiency for a grain size 0.0559$\mu$m with 0.515$\mu$m
ice thickness composed of $\rm H_{2}O$, $\rm NH_{3}$, $\rm CO$, $\rm CO_{2}$,
$\rm CH_{4}$, $\rm CH_{3}OH$ in ratio 0.16:0.16:0.16:0.2:0.16:0.16 with 20% of
porosity. The bottom panel shows the relative error between the interpolation
and explicit calculation of efficiency for the icy grain.

###  3.3 Effective medium vs core-mantle Mie theory

The composition of the icy grains can be thought of as (1) a mixture of ices
and refractory material or (2) a multi-layered structure with refractory
material as the core and an ice mantle on top. The opacity calculations for
these two cases can be carried out using either Mie theory with effective
medium theory or core-mantle Mie theory, respectively. However, the
computational effort required for core-mantle solutions is a factor 4 larger
than the other.

Figure 5: Extinction opacity of bare grains (green) and ice covered dust
grains calculated using the effective medium theory (red) and Core-Mantle Mie
theory (blue) for ice powerlaws $q\\!=\\!0$ (solid lines) and $q\\!=\\!1$
(small dashed lines). The ice composition is $\rm H_{2}O$:$\rm NH_{3}$:$\rm
CH_{4}$:Vacuum = 0.51:0.40:0.04:0.05. The bare grain size distribution,
composition, and volume of ice is the same as in Fig. 2 Figure 6: Comparison
of disk properties of the models. Panels A-F show the vertical visual
extinction, radial visual extinction, dust temperature, gas temperature,
spectral energy distribution, and disk mass averaged opacity of the entire T
Tauri disk [in $\rm cm^{2}/g$(refractory material)] respectively for four disk
models. The four models are 1) bare grains, 2) icy grains with constant
thickness of ice layer ($q\\!=\\!0$), 3) icy grains with constant
ice/refractory volume ratio ($q\\!=\\!1$), 4) icy grains with constant ice
thickness ($q\\!=\\!0$) with grain surface reactions. All the models use
steady state chemistry. Figure 7: Contribution of different grain sizes
(refractory + ice) to the disk mean opacity for bare grain model, constant ice
thickness model, constant ice volume fraction, and constant ice thickness with
grain surface chemistry model. The legend at the top indicates the curves
corresponding to different grain size bins.

It can be argued that in the limiting case $q\\!=\\!1$, the coagulated or
shattered grains do not retain a core-mantle structure. Further, it is
debatable whether icy grains are layered or mixed in dense environments (Noble
et al., 2017). For example, ices could be organized as layers starting with a
layer of polar ices such as $\rm H_{2}O$, $\rm NH_{3}$, and $\rm CH_{3}OH$ and
on top a layer of non-polar ices such as $\rm CO$, N2, and O2, creating an
onion-like structure (Öberg, 2009). Ice covered grains can undergo processing
due to thermal processes, ultraviolet photons, and cosmic ray particles,
called energetically processed ice (Boogert et al., 2015), which can mix the
ices in the mantle and as well cause segregation. Impinging energetic
particles and irradiation with X-rays and UV photons may cause ice desorption
(Palumbo et al., 2006; Andersson & van Dishoeck, 2008), mantle explosion
(Ivlev et al., 2015), or the formation of new ice species (Allamandola et al.,
1988; Muñoz Caro et al., 2002; Palumbo et al., 2006), possibly followed by
molecular ejections (Dulieu et al., 2013) which can lead to mixed ices. The
mixed ice components can undergo processing leading to ice restructuring and
segregation (Fayolle et al., 2011; Öberg et al., 2009a). Further, thermal and
radiative processes can also form larger and more complex molecules (Garrod et
al., 2008; Gerakines et al., 1996; Öberg et al., 2009b; Muñoz Caro et al.,
2002). In the present model, neither ice segregation nor restructuring is
considered. Moreover, using Eq. (6) does not take into account the ice
features that result from interactions between different ice species, for
example $\rm CO$-$\rm H_{2}O$ at 4.68 $\mu$m, 4.647 $\mu$m (Tielens et al.,
1991; Sandford et al., 1988), $\rm NH_{3}$-$\rm H_{2}O$ at 3.4 $\mu$m (Bertie
& Morrison, 1980). Figure 5 shows the ice opacities calculated for the
limiting ice powerlaws with effective medium theory and core-mantle Mie theory
(here, the mantle is a composed of a mixture of ices). The differences between
the opacity values calculated by these two approaches are small. Due to these
reasons our opacity calculations are based on effective medium theory.

##  4 Impact of ice opacities on the disk models

In the following, we discuss four protoplanetary disk models with identical
stellar, disk shape and mass, dust settling, and bare grain parameters, see
Appendix B. The models only differ with respect to the treatment of the ice
opacities and grain surface chemistry. All four models use steady state
chemistry and the fourth model includes grain surface chemistry. The four
models are:

  * •

bare grain: bare grain opacity.

  * •

const ice thickns: ice opacities with a constant ice thickness, $q\\!=\\!0$
(Eq. 3).

  * •

const ice vol frac: ice opacities with a constant ice volume fraction,
$q\\!=\\!1$ (grain size dependent ice thickness).

  * •

const ice thicknss + surf chem: same as the constant ice thickness model, that
is, $q\\!=\\!0$, but includes grain surface chemistry.

The position-dependent icy grain opacity calculations are based on the
abundances of ice species as provided by the chemical network at each point in
the disk. The iteration between chemistry, opacity calculation, and continuum
radiative transfer, as described in Sect. 2, is continued until the opacity
computations converge to a tolerance of less than 5%.

###  4.1 Extinction and opacity

The inclusion of the ice opacities significantly changes the extinction
properties of the disk in the model, particulary in the midplane, where ices
are abundant. Fig. 6 shows the vertical visual extinction, radial visual
extinction, dust temperature, gas temperature, spectral energy distribution,
and disk mass averaged opacity of the entire T Tauri disk [in $\rm
cm^{2}/g$(refractory material)] for the four disk models. The assumptions of
constant ice thickness ($q\\!=\\!0$) and constant ice volume fraction
($q\\!=\\!1$) modifies the grain size distribution differently (Fig. 1), which
in turn affects the mean disk opacity. Panel F of Fig. 6 shows that including
ice opacities significantly increases the disk mean opacity compared to the
bare grain opacity (black) model. The constant ice thickness model (red) shows
an increase in the opacity at shorter wavelengths (¡3 $\mu$m) by more than an
order of magnitude. At longer wavelengths (¿100 $\mu$m), the opacity change
with respect to bare grain model is negligible. The mean opacity of the
constant ice volume fraction model (blue) is larger than the bare grain model
at all wavelengths. Particularly, at shorter wavelengths the increase in
opacity is by a factor $\sim$2 and at longer wavelengths by a factor of
$\sim$1.3. At intermediate wavelengths, the opacities of both constant ice
thickness and constant ice volume fraction models are of the same order of
magnitude. However, the constant ice volume fraction model produces weaker ice
features. The model with constant ice thickness and grain surface chemistry
(green) shows that an increase in the ice abundance (due to grain surface
reactions) leads to a larger opacity compared to all the other models in the
mid-IR wavelengths. For a given volume of ice at one location in a disk, the
effect of the ice thickness to grain size relation on the opacity can also be
seen in Fig. 5. The difference in the opacities between the models is an
interplay of changing grain size distribution due to the ice thickness to
grain size relation and the absolute abundance of ice due to chemistry which
is explained in the following paragraph.

The model with constant ice thickness produces ice features that are stronger
than for the model with constant ice volume fraction by an order of magnitude.
However, the ice abundances are similar in these models. This implies that ice
distribution across the grain sizes become important for the ice features.
Figure 7 shows the contribution of different grain sizes to the disk mean
opacity for the four models. Grains of particular sizes ($a$) have opacity
maxima at wavelengths regions corresponding to $2\pi a/\lambda\sim 1$. At
wavelengths corresponding to $2\pi a/\lambda>>1$, the Mie extinction
efficiencies reach a constant value ($Q_{ext}\to 2$). At wavelengths
corresponding to $2\pi a\lambda<<1$, the extinction is dominated by absorption
($Q_{ext}\simeq Q_{abs}$) and falls off inversely proportional to the
wavelength ($Q_{abs}\propto\lambda^{-1}$) (Krugel, 2002). This behavior is
seen in the bare grain model (left). For the assumed refractory grain size
distribution, the opacity in the wavelength range 2 $\mu$m-100 $\mu$m is
dominated by grains of sizes 1.0 $\mu$m$\leq$a$<$10 $\mu$m and 10
$\mu$m$\leq$a$<$100 $\mu$m. Further, the small grains (¡1 $\mu$m) have strong
silicate (and ice for other models) features. In the constant ice thickness
model, small grains (¡1 $\mu$m) have large volume fractions of ice ($\sim$90%)
due to larger surface area. This leads to change in grain size distribution as
shown in Fig. 1, where there is a pile up of grains around the size
corresponding to ice thickness and removal of grains smaller than the ice
thickness. Due to these changes, the opacity contribution of grains smaller
than local ice thickness drastically reduces ($\sim$1 order) and that of
grains of size corresponding to the ice thickness drastically increases
($\sim$1 order) (Fig. 7). That is, the main opacity carriers in the infrared
wavelengths changes from 1.0 $\mu$m$\leq$a$<$100 $\mu$m to 0.1
$\mu$m$\leq$a$<$10 $\mu$m. This effect is more pronounced in the constant ice
thickness model with grain surface chemistry (right) in which the relevant
grain sizes are 0.1 $\mu$m$\leq$a$<$1 $\mu$m. In contrary, in the constant ice
volume fraction model most of the ice resides on larger grains (¿100 $\mu$m)
which do not contribute to opacity in the near- and mid-infrared wavelengths
and the change in size distribution is uniform across grain sizes. Thus,
models with constant ice thickness show stronger ice features. The opacity
curve of the bare grain model shows two silicate features at 10 $\mu$m and 20
$\mu$m. The constant ice volume fraction model shows small features around
2.95 $\mu$m ($\rm NH_{3}$), 3 $\mu$m ($\rm H_{2}O$), 4 $\mu$m ($\rm CO_{2}$),
26 $\mu$m ($\rm NH_{3}$), and the silicate features. Due to the changes in the
size distribution, the models with constant ice thickness (with and without
grain surface chemistry) show significantly enhanced ice features, prominently
showing $\rm H_{2}O$, $\rm NH_{3}$, $\rm CO_{2}$, $\rm CO$, and $\rm CH_{4}$.
Hence, the average opacity of the disk strongly depends on 1) the ice power
law, $q$, in Eq. (3) and 2) total abundance of ices.

Panels A and B show the effect of the ice opacities on the radial and vertical
visual extinctions. When ice opacities are included, the vertical extinction
increases significantly toward the midplane, but only once the disk has become
optically thick. This increase is more than an order of magnitude in the
models that use a constant ice thickness (in which the opacity increases
significantly as well) closer to the midplane, like an opacity wall or
barrier. For example, the vertical visual extinction at 30 au in the midplane
changes from about 10 to 100 when $q\\!=\\!0$ ice opacities are considered.
The observable disk extent can be defined as the vertical optical depth be
equal to unity at optical wavelengths. This observable disk extent increases
from about 120 au for the bare grain model to about 300 au for the models with
constant ice thickness. This behavior is mostly a size effect, because ice
formation increases the sizes of the small grains in particular, which are the
carriers of the opacity at optical wavelengths. However, the impact of the icy
opacities in the $q\\!=\\!1$ model is less pronounced. Hence, the visual
extinction profile in a disk strongly depends on the relation between the ice
thickness and grain size, that is, the ice power law index $q$ in Eq. (3).

Figure 8: Spectral indices for $\rm 1-4\ \mu m$ (left), $\rm 35-300\ \mu m$
(middle) and $\rm 0.85-1.30\ mm$ (right) for the four models presented in
Sect. 4. Different markers represent the spectral indices for different disk
inclinations as indicated in the legend on the right.

###  4.2 Spectral energy distribution and temperature profile

Although the mean opacities in the disk increase drastically when ice coating
is included, the main disk properties such as temperature, pressure, and
chemical abundances, etc., do not change significantly in the models. This is
because the ice only builds up in the (radially and vertically) optically
thick and shielded regions in the disk. In particular, there is almost no ice
in the model at locations where stellar or interstellar UV photons can
penetrate, even if the local dust temperatures are low. This is because of the
UV photodesorption process, which removes ice species as they interact with UV
photons. The changes of opacity in the optically thick regions have no
significant influence on the local mean intensities ($J_{\nu}$), as
$J_{\nu}\approx B_{\nu}(T_{\rm d})$ is valid anyway ($B_{\nu}$ is the Planck
or black-body radiation), and hence there is no significant influence on the
radiation that escapes from these regions. Panels C and D of Fig. 6 show the
computed dust and gas temperatures in the four disk models. The surface layer
of the disk re-emits or scatters radiation received from the star toward the
observer and toward the disk midplane (Natta, 1993). The latter largely
determines the temperatures in the disk (Chiang & Goldreich, 1997; Chiang et
al., 2001). Since the ice is formed in regions that are already optically
thick, there is almost no ice in the upper regions of the disk that can
scatter incoming radiation. Thus, the temperature structure does not change
much when position dependent ice opacities are included. For example, the
largest change in temperature at any grid point compared to the bare grain
model is 18%; this is for the model with grain surface chemistry. The overall
temperature changes are less than 1.5% for the models without grain surface
chemistry and less than 5% for the model with grain surface chemistry.

The spectral energy distributions (SEDs) of the four models are presented in
panel E of Fig. 6 for disk inclinations 0o, 45o, 75o, and 85o. Since the
visual extinction and the albedo in the disk surface, and the temperature in
the disk largely remain unaffected, we do not see substantial change in the
shape of the SED. However, at higher inclinations due to large optical depths
in the midplane, the flux density at shorter wavelengths (¡3$\mu$m) for the
models with constant ice thickness is about an order of magnitude smaller
compared to the bare grain model. Similarly, the increase in opacity at longer
wavelengths for the model with constant ice volume fraction (for example, the
opacity at 10${}^{\rm 3}\,\mu$m is 1.17 times that of bare grain model)
results in slightly lower flux density in the longer wavelength region of the
SED. This is encapsulated in Fig. 8 which shows the spectral indices ($\rm
n_{a-b}$) for $\rm 1-4\ \mu m$, $\rm 35-300\ \mu m$, and $\rm 0.85-1.3\ mm$
for the four models at different disk inclinations. The spectral index is
defined by the following relation:

$n_{a-b}=\frac{log_{10}(\lambda_{a}F_{a})-log_{10}(\lambda_{b}F_{b})}{log_{10}(\lambda_{a})-log_{10}(\lambda_{b})}$
(7)

where $\lambda F$ is the flux density and $\lambda$ is the wavelength. The
three spectral indices probe different regions of the disk. n${}_{\rm
1-4\leavevmode\nobreak\ \mu m}$ corresponds to the hot inner disk surface,
n${}_{\rm 0.85-1.3\leavevmode\nobreak\ mm}$ corresponds to the cooler outer
disk midplane, and n${}_{\rm 35-300\leavevmode\nobreak\ \mu m}$ corresponds to
warm disk regions. The spectral range for the latter is large because multiple
ice features are present between these wavelengths and determination of the
slope becomes difficult for smaller wavelength ranges, see Fig. 9.

The spectral indices of the models at near- to mid-infrared wavelengths
(n${}_{\rm 1-4\leavevmode\nobreak\ \mu m}$) remain constant across the models
for a given inclination. This is because the radiation at these wavelengths is
emitted from the surface layers close to the inner rim of the disk. These
regions are devoid of ices due to thermal desorption. At millimeter
wavelengths, the radiation arises from large part of the disk and hence can be
used to estimate the dust mass, size of the grains, and their evolution in the
disk (Natta et al., 2007; Testi et al., 2014). Thus, spectral indices at
millimeter wavelengths (n${}_{\rm 0.85-1.3\leavevmode\nobreak\ mm}$) are used
to quantify the dust properties (e.g. Woitke et al., 2019; Villenave et al.,
2020). The spectral indices for any given model depends on the inclination
(Fig. 8). In other words, the model with the same dust properties exhibit
different mm-spectral indices at different inclinations. Secondly, for a given
disk inclination, the spectral index at millimeter wavelengths can depend on
the ice opacity and chemistry model. Since the radiation at these wavelengths
arise from the bulk of the outer and cooler disk, the disk mass averaged
opacity can be used to understand the behavior of the spectral indices. Panel
F of Fig. 6 shows that the model with a constant ice volume fraction has a
steeper opacity curve at these wavelengths and the model with a constant ice
thickness has essentially the same opacity curve as the bare grain model. This
is also reflected in the spectral indices at lower inclinations (0o, 45o).
However, at higher inclinations, the orientation effects play a role leading
to certain parts of the disk being less visible and some being more visible.
At mid- to far-infrared wavelengths, the spectral index (n35-300μm) becomes
more negative (steeper slope) for models with ice opacities without grain
surface chemistry. The model with constant ice volume fraction has a lower
spectral index than the model with constant ice thickness at all disk
inclinations. This shows that the ice power law has an effect on the slope of
mid- to far-infrared SED. The model with constant ice thickness with grain
surface chemistry shows the spectral indices similar to that of the model
without grain surface chemistry. The emission region corresponding to these
wavelengths spans across large range of temperatures (due to large wavelength
range), hence it is difficult to qualitatively analyze the complex interaction
of ice opacities and the spectral index. The indices for the model with grain
surface chemistry for disk inclination of 85o are unreliable due to
substantial ice features, hence are not shown. The ice power law determines
the change in grain size distribution (Fig. 1) and the distribution of ice
across the grain sizes. Thus, the ice power law influences the slope of the
SED at longer wavelengths.

Figure 9: Spectral energy distributions and optical depths of ices in
different models. The top panels show the spectral energy distributions of
different models between $\rm 40\ \mu m$ and $\rm 300\ \mu m$. The black
curves indicate a pseudo-continuum used to obtain the optical depths of the
ice features. In each panel, the SED of the different models are arbitrarily
shifted to provide better visualization. The bottom panels show optical depths
of the ice features in the SEDs in the top panels. Panel A shows the ice
features in the four models at a disk inclination of 75o. Panel B shows the
ice features for the model with constant ice thickness without grain surface
chemistry at different disk inclinations. Note that the orange curve in each
panel is the same model. The ice features are marked using colored rectangles
in the optical depth plots. Figure 10: Simulated images and spatially resolved
spectral flux densities of the disk models presented in this work. The top
panels show simulated images (intensity maps) at 26 $\mu$m and inclination of
75o. The cyan rectangles (0.2${}^{\prime\prime}\times$ 0.3′′) indicates the
region whose spectral flux density is plotted in bottom panel. The spectral
flux density plot also shows the absorption coefficients of the ices and
silicates.

###  4.3 Ice emission and absorption features

Figure 11: Optical depth of the 3 $\mu$m $\rm H_{2}O$ ice feature derived from
the spectral flux density obtained from model 2 in Fig. 10. The optical depth
is shown in dark gray thick line and the absorption coefficients of ice
species that have a feature in the plotted wavelength range are shown as
colored thin lines.

The models produce ice absorption and emission features in the SED. At
inclinations lower than the edge-on case, the SEDs do not show any ice
absorption features in any of the four models. However, weak emission features
at far-infrared wavelengths can be observed corresponding to $\rm H_{2}O$,
$\rm NH_{3}$ ice, and $\rm CO_{2}$ ice. Close to edge-on ($\rm 85^{o}$),
models with constant ice thickness show ice absorption features. Particularly,
the model without grain surface chemistry shows only a strong $\rm CO$ ice
absorption feature and a weak $\rm NH_{3}$ ice absorption feature at 4.6
$\mu$m and 26 $\mu$m respectively. The model with grain surface chemistry
shows significantly more ice absorption features (and stronger) in addition to
those seen in previous model, such as $\rm H_{2}O$, $\rm CO_{2}$, and $\rm
CH_{3}OH$. This model also shows stronger ice emission features at far-IR
wavelengths. This shows that the chemical network strongly affects the
strength of the ice absorption and emission features. Therefore, observations
of such ice features can lay strong constraints on the type of chemistry in
disks.

Figure 9 shows the zoomed-in SED (top panels) and optical depths (bottom
panels) of the ice emission features (at $\rm 75^{o}$ inclination) from $\rm
40-300\ \mu m$ for the models discussed before (Sect. 4). The black curves in
the SED plots in all the top panels show the pseudo continuum used to obtain
the optical depths of the ice features. This continuum is obtained by spline
interpolation using python library scipy (Virtanen et al., 2020) and then
correcting spline of other models for systematic error produced by the spline
for bare grain model. The top and bottom plots in Panel A show the SED and
optical depth of the four models. Due to higher ice abundances, the model with
grain surface chemistry produces stronger features (by about a factor 2 in
optical depth) than the model without grain surface chemistry. The peaks of
the optical depths of emission features for the constant ice thickness and the
constant ice volume fraction models are quite similar, however, the constant
ice volume fraction model shows larger optical depths (by a factor less than
2) between the peaks of these features. The figure also shows the ices that
correspond to the features, the prominent ices being $\rm H_{2}O$, $\rm
NH_{3}$, $\rm CO_{2}$, and $\rm CH_{3}OH$. Unlike other ices, $\rm CH_{3}OH$
ice has a broad feature that spans from about 50 $\mu$m to 160 $\mu$m and
hence does not show up as a single peak instead contributes to the optical
depth across the wavelength range. Panel B of Fig. 9 shows that the optical
depths are quite similar across the inclinations except for the close to edge-
on case, where some of the features (such as $\rm 47\leavevmode\nobreak\ \mu
m$, $\rm 70\leavevmode\nobreak\ \mu m$ $\rm H_{2}O$ and $\rm NH_{3}$ ice
features) can be seen in absorption. Thus, the strength of the ice features
depend on the chemistry producing the ices, the relation governing ice
distribution over different grain sizes and the disk inclination.

The SED shows absorption features only at high disk inclinations. Hence, it is
interesting to explore the observability of absorption features when the disk
is spatially resolved. Figure 10 shows images of the four models at 26 $\mu$m
in the top panels. The bottom of the figure shows the spectral flux density of
these models obtained from the region marked by the cyan rectangles. The cyan
rectangles are $\rm 0.2^{\prime\prime}\times 0.3^{\prime\prime}$ in
size222which can be similar to binning 6 pixels of JWST NIRSpec IFU pixels
(which are $\rm 0.1^{\prime\prime}\times 0.1^{\prime\prime}$ each). This is a
simple study to estimate the fluxes for spatially resolved observations and is
not an attempt to simulate any specific instrument. Hence, point spread
function, noise, errors, sensitivity, and saturation limits are not
considered. Simulating observations and identifying best observation
strategies are left to a future work.

The images in the top panels show that inclusion of ice opacities produces a
distinct band of absorption around the midplane, subsequently called the dark
lane. This dark lane is more prominent for models with constant ice thickness
and also depends on the chemical network. These differences reflect the
increase in the extinction in the disks as discussed in Sect. 4. Panel A of
Fig. 6 shows that the inclusion of ice opacities increases the extinction
radially, more strongly for models with constant ice thickness compared to
model with constant ice volume fraction. The ice power law has a strong effect
on the appearance of the disks, particularly on the size of the dark lane. The
spectral flux density (in the bottom panel) extracted in the dark lane shows
that ice absorption features can be seen by spatially resolving the disks at
inclinations less than the edge-on case. Further, the constant ice volume
fraction model shows substantially weaker features compared to the models that
consider a constant ice thickness. The absorption coefficients of ices and the
refractory silicate used are also shown at the bottom of the figure to help
identify the features observed in the spectra. More information regarding the
optical constants used in this paper is described in Appendix A. The models
show features corresponding to five of the six ices included in the opacity
calculation. Though the local abundance of $\rm CH_{4}$ ice can be as large as
$1.4\cdot 10^{-4}$ in the disk, no feature corresponding to it is observed.
This is because, though the local $\rm CH_{4}$ abundance can be large, the
column density along the line of sight of the spectra extracted in Fig. 10 for
other ices such as $\rm H_{2}O$, $\rm NH_{3}$ and $\rm CH_{3}OH$ are larger
due to the spatial extent of these ices (see appendix). Moreover, the $\rm
CH_{4}$ features lie in the wavelength regions with features from more
dominating ices. This can also be seen in the disk mean opacity of the model
with grain surface chemistry as shown in panel F of Fig. 6. The model with
grain surface chemistry shows stronger absorption features than the model
without grain surface chemistry by about 2 orders of magnitude. Ice emission
features are also present at longer wavelengths but are substantially weaker
(by several orders of magnitude) compared to the absorption features. Hence,
spatially resolving the disk produces significant ice absorption and emission
features that are not seen in the disk integrated spectra, particularly strong
in the absorption by few orders of magnitude.

Figure 11 shows the optical depth of the $\rm 3\ \mu m$ ice feature obtained
from spatially resolved spectral density plot of the model with constant ice
thickness from Fig. 10, along with the absorption coefficients of different
ice species that have a feature in that wavelength range. The 3 $\mu$m profile
is a combination of $\rm H_{2}O$, $\rm NH_{3}$, and $\rm CH_{3}OH$ ices. The
sharp peak at 2.95 $\mu$m in the profile is due to $\rm NH_{3}$ ice. $\rm
H_{2}O$ ice dominates the profile of the 3 $\mu$m feature, however, the longer
wavelength side (¿3.3 $\mu$m) shows a contribution of $\rm CH_{3}OH$ ice.
NIRSpec IFU should be able to spatially resolve the disks and hence we can
expect observations of such features from JWST.

##  5 Discussion

Ice accretion on the grains modifies the grain size distribution. In Sect. 3.1
we showed that a constant ice thickness leads to a deviation of the size
distribution from the powerlaw and a constant ice volume fraction increases
the minimum and maximum grain sizes. Further, these changes can be enhanced by
models that produce more ice, such as the model with grain surface chemistry.
The emission at (sub) millimeter wavelengths probe the bulk of the dust in the
disk. This corresponds to the large dust grains which are settled in the disk
midplane. In these environments grain evolution processes such as grain growth
are important as they affect the emission at these wavelengths. Natta et al.
(2007) show that the millimeter slope is strongly affected by the maximum
grain size and the dust distribution powerlaw. Particularly, when the maximum
grain size is of the order of millimeters, the slope is very sensitive to the
changes in maximum grain size. Since, the emission in these wavelengths can be
used to quantify the dust properties and their evolution, it is important to
consider the effect of ice on the grain size distribution and the opacity
while interpreting these spectral indices.

Ices have been detected in protoplanetary disks, as discussed in Sect. 1.
These detections are from the total disk integrated fluxes (SED) and are not
spatially resolved. Our models show weak emission features at far-IR
wavelengths at all inclinations, which can be, however, very weak and barely
observable for some model setups. Absorption features can only be seen in the
edge-on orientation in the total disk spectra. Even when seen edge-on, most of
the absorption features are seen only when grain surface reactions are
included, for example, the $\rm 3\ \mu m$ $\rm H_{2}O$ feature. However,
spatially resolving the disk produces a multitude of absorption features
across the spectral range.

In contrary, the $\rm 3\ \mu m$ water ice feature has been reported in total
disk flux of some objects (mostly edge-on) such as HK Tau B (Terada et al.,
2007), d132-1832, d216-0939 (Terada et al., 2012). Further, emission features
are observed in a few protoplanetary disks (one such example is VW Cha,
McClure et al., 2015). The presence of emission features and absence of
absorption features in our models can be attributed to the limited parameters
explored. Objects differ, for example in terms of the parent star, disk mass,
the inclination, the flaring angle, inner and outer disk radii, gaps,
silhouette disks, etc. All these parameters can strongly affect the ice
features. Consequently, each disk has to be modeled separately and our models
are not representative of every disk. Further, dynamic physical processes such
as vertical and radial transport are not included in our models and can
significantly alter the ice features. For example, vertical mixing can
significantly change the ice composition in the disk (Furuya & Aikawa, 2014)
and readily produces ice features that are stronger compared those from our
models by factors of a few in terms of optical depth (Woitke et al. in prep).
Similarly, models with gaps can also show strong ice features in the SED. The
goal of this paper is to present the consistent position dependent opacities,
and exploring such parameters and simulating observations is left to future
work. For example, the disks with ice features in the far-IR reported by
McClure et al. (2015) can be modeled with our ice opacity approach.

Ballering et al. (2021) also modeled ices in disks and they find absorption
features of multiple ices in the total disk flux even at inclinations lower
than the edge-on case. However, these models differ from our models
significantly and hence a direct comparison is not possible. For example,
Ballering et al. (2021) do not consider the interstellar UV field and use a
cosmic ray ionization rate that is an order of magnitude lower than our
models. These assumptions can significantly affect the ice features. The
vertical extent of ice in the disk is largely determined by the interstellar
UV field, because this radiation can photodesorb the ices back into the gas
phase. This in turn inhibits the presence of ice in the emitting regions of
the disk. The ice destruction mechanism by photodesorption is consistently
treated in our models, with chemistry-dependent ice opacities and fully
coupled UV radiative transfer. Further, Ballering et al. (2021) use two dust
populations with significantly different settling conditions and the treatment
of ice opacities is not fully consistent with local conditions in each grid
point of the model.

The strengths of the ice features strongly depend on the ice power law and the
chemistry network. Spatially resolving the disk presents multiple ice features
and can probe the vertical and radial distribution of ice in the disk. In the
end, future observations will allow us to test these various models and to
refine the assumptions.

##  6 Conclusions

A computationally efficient way of computing position-dependent ice and bare
grain opacities in protoplanetary disks has been implemented in PRODIMO, which
consistently couples chemistry, ice formation, opacity computations, and
radiative transfer in protoplanetary disks. The conclusions drawn are the
following:

  1. 1.

A weighted sum of the opacities of grains coated with the different pure ices
and log interpolation of opacities between pre-calculated points in ice
thickness space provides a computationally efficient and sufficiently accurate
method to calculate the icy opacities in protoplanetary disks.

  2. 2.

Locally, the continuum opacity can change significantly by ice formation. The
opacity changes are particularly large at shorter wavelengths, due to the
increased size of the small grains, when assuming a constant ice thickness on
all grains. Opacity changes are less pronounced, and more uniform in
wavelength space, when the grains are assumed to have a constant
ice/refractory volume ratio, which is expected if coagulation and shattering
are important. In both cases, the primary reason for those opacity changes is
the altered size distribution of the icy grains.

  3. 3.

The position-dependent ice opacities do not alter the physical disk structure
significantly, because the ice opacities are significant only in the optically
thick disk regions ($A_{V}\\!\ga\\!10$) in our models.

  4. 4.

The spectral energy distributions (SEDs) of the disks show faint ice emission
features in the far-infrared regions beyond 35 $\mu$m for all our models with
ice opacities. However, absorption features are seen only at high
inclinations, close to edge-on conditions. In particular, the 3 $\mu$m water
ice feature is seen only in case of the model with grain surface chemistry in
close to edge-on orientation.

  5. 5.

The assumption made on how the ice is distributed across the grain size
distribution (ice power law) influences the far-infrared and millimeter slope
of the SED.

  6. 6.

The strengths of the ice features are influenced by the ice power law (about a
factor of 2 in optical depths). The number of ice features and their strength
increase (a factor 2 in optical depth) in models having grain surface
chemistry, which is directly related to an increase in ice abundance.

  7. 7.

Our models show a prominent dark lane in resolved images. This dark lane is
caused by icy grains situated just above the disk midplane, and the size of
the dark lane depends on the assumptions used to calculate the modified dust
size distribution and the chemical network.

  8. 8.

To detect the ice absorption features, extracting spatially resolved spectra
from the dark lane is more promising.

The models presented here are useful to model existing (e.g. Herschel) and
upcoming observations (e.g. JWST) to better understand the ice composition in
disks. And further improvements, such as transport processes, will help us
understand the role of these processes in disk environments.

###### Acknowledgements.

CHR acknowledge the support of the the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) - 325594231. CHR is grateful for support from the
Max Planck Society.

## References

  * Allamandola et al. (1999) Allamandola, L. J., Bernstein, M. P., Sandford, S. A., & Walker, R. L. 1999, Space Sci. Rev., 90, 219
  * Allamandola et al. (1988) Allamandola, L. J., Sandford, S. A., & Valero, G. J. 1988, Icarus, 76, 225
  * Andersson & van Dishoeck (2008) Andersson, S. & van Dishoeck, E. F. 2008, A&A, 491, 907
  * Ballering et al. (2021) Ballering, N. P., Cleeves, L. I., & Anderson, D. E. 2021, ApJ, 920, 115
  * Bergin et al. (2007) Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck, E. F. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 751
  * Bertie & Morrison (1980) Bertie, J. E. & Morrison, M. M. 1980, J. Chem. Phys., 73, 4832
  * Boogert & Ehrenfreund (2004) Boogert, A. C. A. & Ehrenfreund, P. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 309, Astrophysics of Dust, ed. A. N. Witt, G. C. Clayton, & B. T. Draine, 547
  * Boogert et al. (2015) Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541
  * Boogert et al. (2008) Boogert, A. C. A., Pontoppidan, K. M., Knez, C., et al. 2008, ApJ, 678, 985
  * Brooke et al. (1999) Brooke, T. Y., Sellgren, K., & Geballe, T. R. 1999, ApJ, 517, 883
  * Brown & Bolina (2007) Brown, W. A. & Bolina, A. S. 2007, MNRAS, 374, 1006
  * Bruderer et al. (2012) Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91
  * Chiang & Goldreich (1997) Chiang, E. I. & Goldreich, P. 1997, ApJ, 490, 368
  * Chiang et al. (2001) Chiang, E. I., Joung, M. K., Creech-Eakman, M. J., et al. 2001, ApJ, 547, 1077
  * Chiar et al. (2002) Chiar, J. E., Adamson, A. J., Pendleton, Y. J., et al. 2002, ApJ, 570, 198
  * Chiar et al. (2000) Chiar, J. E., Tielens, A. G. G. M., Whittet, D. C. B., et al. 2000, ApJ, 537, 749
  * Collings et al. (2004) Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133
  * D’Alessio et al. (2006) D’Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., & Servín, H. 2006, ApJ, 638, 314
  * Dartois & d’Hendecourt (2001) Dartois, E. & d’Hendecourt, L. 2001, A&A, 365, 144
  * Dartois et al. (2002) Dartois, E., d’Hendecourt, L., Thi, W., Pontoppidan, K. M., & van Dishoeck, E. F. 2002, A&A, 394, 1057
  * Dorschner et al. (1995) Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503
  * Dulieu et al. (2013) Dulieu, F., Congiu, E., Noble, J., et al. 2013, Scientific Reports, 3, 1338
  * Edridge (2010) Edridge, J. L. 2010, PhD thesis, University of London, University College London (United Kingdom)
  * Fayolle et al. (2011) Fayolle, E. C., Öberg, K. I., Cuppen, H. M., Visser, R., & Linnartz, H. 2011, A&A, 529, A74
  * Furuya & Aikawa (2014) Furuya, K. & Aikawa, Y. 2014, ApJ, 790, 97
  * Garrod et al. (2008) Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, ApJ, 682, 283
  * Gerakines et al. (1996) Gerakines, P. A., Schutte, W. A., & Ehrenfreund, P. 1996, A&A, 312, 289
  * Gibb et al. (2004) Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A. G. G. M. 2004, ApJS, 151, 35
  * Gibb et al. (2000) Gibb, E. L., Whittet, D. C. B., Schutte, W. A., et al. 2000, ApJ, 536, 347
  * Goto et al. (2018) Goto, M., Bailey, J. D., Hocuk, S., et al. 2018, A&A, 610, A9
  * Herrero et al. (2010) Herrero, V. J., Gálvez, Ó., Maté, B., & Escribano, R. 2010, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 12, 3164
  * Honda et al. (2009) Honda, M., Inoue, A. K., Fukagawa, M., et al. 2009, ApJ, 690, L110
  * Honda et al. (2016) Honda, M., Kudo, T., Takatsuki, S., et al. 2016, ApJ, 821, 2
  * Hudgins et al. (1993) Hudgins, D. M., Sandford, S. A., Allamandola, L. J., & Tielens, A. G. G. M. 1993, ApJS, 86, 713
  * Ivlev et al. (2015) Ivlev, A. V., Röcker, T. B., Vasyunin, A., & Caselli, P. 2015, ApJ, 805, 59
  * Kamp et al. (2018) Kamp, I., Scheepstra, A., Min, M., Klarmann, L., & Riviere-Marichalar, P. 2018, A&A, 617, A1
  * Kamp et al. (2017) Kamp, I., Thi, W. F., Woitke, P., et al. 2017, A&A, 607, A41
  * Kamp et al. (2010) Kamp, I., Tilling, I., Woitke, P., Thi, W. F., & Hogerheijde, M. 2010, A&A, 510, A18
  * Krugel (2002) Krugel, E. 2002, The physics of interstellar dust (CRC Press)
  * Kuroiwa & Sirono (2011) Kuroiwa, T. & Sirono, S.-i. 2011, ApJ, 739, 18
  * Martonchik et al. (1984) Martonchik, J. V., Orton, G. S., & Appleby, J. F. 1984, Appl. Opt., 23, 541
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * McClure et al. (2015) McClure, M. K., Espaillat, C., Calvet, N., et al. 2015, ApJ, 799, 162
  * McClure et al. (2012) McClure, M. K., Manoj, P., Calvet, N., et al. 2012, ApJ, 759, L10
  * Merrill et al. (1976) Merrill, K. M., Russell, R. W., & Soifer, B. T. 1976, ApJ, 207, 763
  * Min et al. (2016a) Min, M., Bouwman, J., Dominik, C., et al. 2016a, A&A, 593, A11
  * Min et al. (2005) Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909
  * Min et al. (2016b) Min, M., Rab, C., Woitke, P., Dominik, C., & Ménard, F. 2016b, A&A, 585, A13
  * Muñoz Caro et al. (2002) Muñoz Caro, G. M., Meierhenrich, U. J., Schutte, W. A., et al. 2002, Nature, 416, 403
  * Natta (1993) Natta, A. 1993, ApJ, 412, 761
  * Natta et al. (2007) Natta, A., Testi, L., Calvet, N., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 767
  * Noble et al. (2013) Noble, J. A., Fraser, H. J., Aikawa, Y., Pontoppidan, K. M., & Sakon, I. 2013, ApJ, 775, 85
  * Noble et al. (2017) Noble, J. A., Fraser, H. J., Pontoppidan, K. M., & Craigon, A. M. 2017, MNRAS, 467, 4753
  * Nomura & Millar (2005) Nomura, H. & Millar, T. J. 2005, A&A, 438, 923
  * Öberg (2009) Öberg, K. I. 2009, PhD thesis, Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands
  * Öberg & Bergin (2021) Öberg, K. I. & Bergin, E. A. 2021, Phys. Rep, 893, 1
  * Öberg et al. (2011) Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2011, ApJ, 740, 109
  * Öberg et al. (2009a) Öberg, K. I., Fayolle, E. C., Cuppen, H. M., van Dishoeck, E. F., & Linnartz, H. 2009a, A&A, 505, 183
  * Öberg et al. (2009b) Öberg, K. I., Garrod, R. T., van Dishoeck, E. F., & Linnartz, H. 2009b, A&A, 504, 891
  * Ossenkopf & Henning (1994) Ossenkopf, V. & Henning, T. 1994, A&A, 291, 943
  * Palumbo et al. (2006) Palumbo, M. E., Baratta, G. A., & Spinella, F. 2006, Memorie della Societa Astronomica Italiana Supplementi, 9, 192
  * Palumbo & Strazzulla (1993) Palumbo, M. E. & Strazzulla, G. 1993, A&A, 269, 568
  * Pontoppidan et al. (2008) Pontoppidan, K. M., Boogert, A. C. A., Fraser, H. J., et al. 2008, ApJ, 678, 1005
  * Pontoppidan et al. (2005) Pontoppidan, K. M., Dullemond, C. P., van Dishoeck, E. F., et al. 2005, ApJ, 622, 463
  * Rab et al. (2018) Rab, C., Güdel, M., Woitke, P., et al. 2018, A&A, 609, A91
  * Rocha & Pilling (2014) Rocha, W. & Pilling, S. 2014, Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy, 123, 436
  * Sandford & Allamandola (1988) Sandford, S. A. & Allamandola, L. J. 1988, Icarus, 76, 201
  * Sandford et al. (1988) Sandford, S. A., Allamandola, L. J., Tielens, A. G. G. M., & Valero, G. J. 1988, ApJ, 329, 498
  * Semenov & Wiebe (2011) Semenov, D. & Wiebe, D. 2011, ApJS, 196, 25
  * Smith et al. (1989) Smith, R. G., Sellgren, K., & Tokunaga, A. T. 1989, ApJ, 344, 413
  * Terada & Tokunaga (2012) Terada, H. & Tokunaga, A. T. 2012, ApJ, 753, 19
  * Terada & Tokunaga (2017) Terada, H. & Tokunaga, A. T. 2017, ApJ, 834, 115
  * Terada et al. (2007) Terada, H., Tokunaga, A. T., Kobayashi, N., et al. 2007, ApJ, 667, 303
  * Terada et al. (2012) Terada, H., Tokunaga, A. T., Pyo, T.-S., et al. 2012, AJ, 144, 175
  * Testi et al. (2014) Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339
  * Thi et al. (2020) Thi, W. F., Hocuk, S., Kamp, I., et al. 2020, A&A, 634, A42
  * Thi et al. (2011) Thi, W. F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711
  * Tielens et al. (1991) Tielens, A. G. G. M., Tokunaga, A. T., Geballe, T. R., & Baas, F. 1991, ApJ, 381, 181
  * Villenave et al. (2020) Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Warren (1984) Warren, S. G. 1984, Appl. Opt., 23, 1206
  * Warren (1986) Warren, S. G. 1986, Appl. Opt., 25, 2650
  * Whittet et al. (1983) Whittet, D. C. B., Bode, M. F., Longmore, A. J., Baines, D. W. T., & Evans, A. 1983, Nature, 303, 218
  * Whittet et al. (1998) Whittet, D. C. B., Gerakines, P. A., Tielens, A. G. G. M., et al. 1998, ApJ, 498, L159
  * Woitke et al. (2019) Woitke, P., Kamp, I., Antonellini, S., et al. 2019, PASP, 131, 064301
  * Woitke et al. (2009) Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103
  * Zubko et al. (1996) Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321

##  Appendix A Ice optical constants

The source of optical data, the measurement temperature, and structure of six
types of ice species and refractory material used are tabulated in Table 1.
Ice structures change with temperature, so does the refractive index and
absorption coefficient. Though the abundances of the different ices follow the
local thermo-chemical state of the disk, the optical constants used are not
consistent with temperature. Table 1 summarizes the state and temperature of
ices at which the optical data is obtained. It should be noted that the ice
optical data do not cover the entire spectrum (0.0912 $\mu$m - 104 $\mu$m).
Hence, the data has to be expanded extensively for all species to enable
studying all the ice features. For absence of optical data corresponding to a
wavelength $\lambda$ which is bounded by optical data at $\lambda_{1}$ and
$\lambda_{2}$ (such that $\lambda_{1}<\lambda<\lambda_{2}$), we perform a log
interpolation to estimate the optical data at $\lambda$. In an unbounded case
(extrapolation), for longer wavelength side ($\lambda>\lambda_{1}$) we assume:
$n(\lambda)=n(\lambda_{1})$ and
$k(\lambda)=k(\lambda_{1})\cdot\lambda_{1}/\lambda$; for shorter wavelength
side ($\lambda<\lambda_{2}$) we assume: $n(\lambda)=n(\lambda_{2})$ and
$k(\lambda)=k(\lambda_{2})$. The optical constants used are shown in Fig. 12
and Fig. 13.

Figure 12: Refractive indices (n) of ices species used in the models. Figure 13: Absorption coefficients (k) of ices species used in the models. Table 1: Optical data of different refractory and ice species. Ice Species |  Temperature (K) |  Type |  Reference
---|---|---|---
$\rm H_{2}O$ |  266 |  Crystalline |  Warren (1984)
$\rm NH_{3}$ |  77-88 |  Crystalline |  Martonchik et al. (1984)
$\rm CO_{2}$ |  Multiple |  Crystalline |  Warren (1986)
$\rm CO$ |  10 |  Amorphous |  Hudgins et al. (1993); Rocha & Pilling (2014)
$\rm CH_{3}OH$ |  10 |  Amorphous |  Hudgins et al. (1993)
$\rm CH_{4}$ |  10 |  Amorphous |  Hudgins et al. (1993)
Carbon |  |  Amorphous |  Zubko et al. (1996)
Silicate (Mg0.7Fe0.3SiO3) |  |  Amorphous |  Dorschner et al. (1995)

##  Appendix B Disk model

We present a fiducial T Tauri disk ProDiMo model. Table 2 presents the disk
properties considered for all the models presented in this paper. For more
details regarding the parameters, meanings and their values, refer Woitke et
al. (2009) and Woitke et al. (2016).

Figure 14: Distribution of the six ice species ($\rm H_{2}O$, $\rm CO$, $\rm
CO_{2}$, $\rm NH_{3}$, $\rm CH_{3}OH$, and $\rm CH_{4}$) in the presented four
models with different ice opacity ($q$) and chemistry settings (See Sect. 4).
Each row corresponds to the different ice species and each column corresponds
to the one of the four models. The dust temperature and vertical visual
extinction contours are also shown with solid and dashed lines respectively.
Table 2: Common disk parameters for the models.

Property |  Symbol |  Value
---|---|---
Stellar Mass |  M∗ |  0.7 M⊙
Effective temperature |  T∗ |  4000 K
Stellar luminosity |  L∗ |  1 L⊙
UV excess |  fUV |  0.01
UV powerlaw index |  pUV |  1.3
X-ray luminosity |  Lx |  1030 erg/s
X-ray emission temperature |  Tx |  2$\times$107 K
Strength of interstellar UV |  $\chi^{\rm ISM}$ |  1
Strength of interstellar IR |  $\chi^{\rm ISM}_{\rm IR}$ |  0
Cosmic ray H2 ionization rate |  $\zeta_{\rm CR}$ |  1.7$\times$10-17 s-1
Disk Mass |  Mdisk |  0.01 M⊙
Dust/gas mass ratio |  $\delta$ |  0.01
Inner disk radius |  Rin |  0.07 au
Tapering-off radius |  Rtap |  100 au
Column density power index |  $\epsilon$ |  1
Reference scale height |  Hg(100 au) |  10 au
Flaring power index |  $\beta$ |  1.15
Minimum dust particle radiusa |  amin |  0.05 $\mu$m
Maximum dust particle radiusa |  amax |  3000 $\mu$m
Dust size dist. power indexa |  apow |  3.5
Turbulent mixing parameter |  $\alpha_{\rm settle}$ |  0.01
Refractory dust composition |  Mg0.7Fe0.3SiO3 |  60%
|  amorph. C |  15%
|  porosity |  25%
PAH abundance rel. to ISM |  fPAH |  0.01
Chemical heating efficiency |  $\gamma^{\rm chem}$ |  0.2
Distance to the observer |  d |  140 pc

(a) Parameters related to refractory dust grains.

Table 3: Adsorption energies of different ice species. Ice species |  Adsorption energy (K) |  Reference
---|---|---
$\rm H_{2}O$ |  4800.0 |  Brown & Bolina (2007)
$\rm NH_{3}$ |  5534.0 |  Collings et al. (2004)
$\rm CO_{2}$ |  2990.0 |  Edridge (2010)
$\rm CO$ |  1150.0 |  Collings et al. (2004)
$\rm CH_{3}OH$ |  4930.0 |  Brown & Bolina (2007)
$\rm CH_{4}$ |  1090.0 |  Herrero et al. (2010)
Figure 15: Mean disk abundance of ice species in the four models. Figure 16:
Ice thickness distribution in the constant ice thickness model without grain
surface chemistry. The green diamond indicates the location of largest ice
thickness (0.61 $\mu$m) Figure 17: Ice thickness distribution in the constant
ice thickness model with grain surface chemistry. The green diamond indicates
the location of largest ice thickness (0.59 $\mu$m) Figure 18: Ice volume
fraction distribution in the constant ice volume fraction model. The green
diamond indicates the location of largest ice volume fraction (0.99)

We use the large DIANA chemical network composed of 235 species with 13
elements. The network is described in detail in Kamp et al. (2017). The
adsorption energies of different ices considered in this paper are given in
Table 3. The grain surface chemistry used in our models is explained in (Thi
et al. 2020), we add two species on top of the large DIANA network: $\rm H\\#$
and $\rm H_{2}\\#$, allowing for a minimal surface chemistry, e.g,
hydrogenation on surfaces. These models only use physisorption and no
chemisorption for the surface chemistry. The purpose of including the surface
chemistry is to assess the impact of different chemical networks, and is not
to explain the ice observations and the compositions in these objects. We also
ran time-dependent chemistry model with grain surface chemistry where the
initial composition was determined by a molecular cloud model. We evolved the
model up to 2Myr. We noticed that the ice abundances, hence the opacities and
the ice features, are very similar to that of the steady state chemistry with
grain surface reactions, hence the model is not included in Sect. 4.

Figure 14 shows the abundances of $\rm H_{2}O$, $\rm NH_{3}$, $\rm CO$, $\rm
CO_{2}$, $\rm CH_{3}OH$, and $\rm CH_{4}$ ices in the four models discussed in
Sect. 4. The figure also shows the dust temperatures and vertical extinction
in the disk models. Each column of panels in the figure shows the abundance of
the ice species in each model. Models with ice opacities show an increase in
water ice in the midplane which extends radially. The outer radii of water ice
region extends from 7 au in bare grain model to 25 au, 10 au, and 400 au in
models with constant ice thickness, constant ice volume fraction, and constant
ice thickness with grain surface chemistry respectively. Fig. 15 shows the
mean abundances of the ices in the models presented in this paper. The models
with ice opacities are generally more abundant in ices, with exception of $\rm
CO$ and $\rm CH_{4}$ ices. The ice abundances are slightly higher in the model
with constant ice thickness than the model with constant ice volume fraction
by about a factor 1.5. The model with grain surface chemistry has several
orders of magnitude higher mean ice abundances, except for $\rm CO$ ice which
is about two orders lower. In general, including ice opacities results in
higher ice abundances except for $\rm CO$ ice. It should be noted that these
are mean abundances of the entire disk, hence, essentially are model
properties and should not be mistaken for the ice abundances that can be
obtained from observations. Figures 16 and 17 show the distribution of ice
thickness in the disk for the models with constant ice thickness without and
with grain surface chemistry respectively. The mean ice thickness for entire
the disk (which includes regions without any ice) in these models are 0.045
$\mu$m and 0.094 $\mu$m respectively. Figure 18 shows the distribution of ice
volume fraction across the disk for the constant ice volume fraction model.
The mean ice volume fraction for the entire disk in this model is 24.2%.
Unlike constant ice thickness models, the constant ice volume fraction model
shows lower ice volume fraction in the midplane and larger ice volume fraction
in the upper midplane. Dust settling leads to an increase in the large grain
population in the midplane. Since ice largely resides on large grains in this
model, the regions with little large grains yet having significant ice
abundances have larger ice volume fractions.