Context. The formation of planetesimals in protoplanetary disks is not wellunderstood. Streaming instability is a promising mechanism to directly form planetesimals from pebblesized particles provided a high enough solidstogas ratio. However local enhancements of the solidstogas ratio are difficult to realize in a smooth disk which motivates the consideration of special disk locations such as the snowline  the radial distance from the star beyond which water can condense into solid ice. BR  Aims In this article we investigate the viability of planetesimal formation by streaming instability near the snowline due to water diffusion and condensation. We aim to identify under what disk conditions streaming instability can be triggered near the snowline. BR  Methods To this end we adopt a viscous disk model and numerically solve the transport equations for vapor and solids on a cylindrical 1D grid. We take into account radial drift of solids gas accretion on to the central star and turbulent diffusion. We study the importance of the backreaction of solids on the gas and of the radial variation of the mean molecular weight of the gas. Different designs for the structure of pebbles are investigated varying in the number and size of silicate grains. We also introduce a semianalytical model that we employ to obtain results for different disk model parameters. BR  Results We find that water diffusion and condensation can locally enhance the ice surface density by a factor 35 outside the snowline. Assuming that icy pebbles contain many micronsized silicate grains that are released during evaporation the enhancement is increased by another factor 2. In this manyseeds model the solidstogas ratio interior to the snowline is enhanced as well but not as much as just outside the snowline. In the context of a viscous disk the diffusioncondensation mechanism is most effective for high values of the turbulence parameter  10SUP3SUP10SUP2SUP. Therefore assuming young disks are more vigorously turbulent than older disks planetesimals near the snowline can form in an early stage of the disk. In highly turbulent disks tens of Earth masses can be stored in an annulus outside the snowline which can be identified with recent ALMA observations.
['accretion', 'accretion disks', 'turbulence', 'methods: numerical', 'planets and satellites: formation', 'protoplanetary disks', 'Astrophysics - Earth and Planetary Astrophysics']
Planetesimal formation near the snowline in or out
{'Planetesimal formation near the snowline: in or out?': 'Djoeke Schoonenberg and Chris W. Ormel  \nAnton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090 GE Amsterdam, The Netherlands e-mail: [d.schoonenberg@uva.nl,c.w.ormel@uva.nl]  \nReceived 4 November 2016 / Accepted 4 February 2017', 'ABSTRACT': "The formation of planetesimals in protoplanetary disks is not well-understood. Streaming instability is a promising mechanism to directly form planetesimals from pebble-sized particles, provided a high enough solids-to-gas ratio. However, local enhancements of the solids-to-gas ratio are di GLYPH<14> cult to realize in a smooth disk, which motivates the consideration of special disk locations such as the snowline - the radial distance from the star beyond which water can condense into solid ice.  \nIn this article we investigate the viability of planetesimal formation by streaming instability near the snowline due to water di GLYPH<11> usion and condensation. We aim to identify under what disk conditions streaming instability can be triggered near the snowline.  \nTo this end, we adopt a viscous disk model, and numerically solve the transport equations for vapor and solids on a cylindrical, 1D grid. We take into account radial drift of solids, gas accretion on to the central star, and turbulent di GLYPH<11> usion. We study the importance of the back-reaction of solids on the gas and of the radial variation of the mean molecular weight of the gas. Di GLYPH<11> erent designs for the structure of pebbles are investigated, varying in the number and size of silicate grains. We also introduce a semi-analytical model that we employ to obtain results for di GLYPH<11> erent disk model parameters.  \nWe find that water di GLYPH<11> usion and condensation can locally enhance the ice surface density by a factor 3-5 outside the snowline. Assuming that icy pebbles contain many micron-sized silicate grains that are released during evaporation, the enhancement is increased by another factor GLYPH<24> 2. In this 'many-seeds' model, the solids-to-gas ratio interior to the snowline is enhanced as well, but not as much as just outside the snowline. In the context of a viscous disk, the di GLYPH<11> usion-condensation mechanism is most e GLYPH<11> ective for high values of the turbulence parameter GLYPH<11> (10 GLYPH<0> 3 -10 GLYPH<0> 2 ). Therefore, assuming young disks are more vigorously turbulent than older disks, planetesimals near the snowline can form in an early stage of the disk. In highly turbulent disks, tens of Earth masses can be stored in an annulus outside the snowline, which can be identified with recent ALMA observations.  \nKey words. accretion, accretion disks - turbulence - methods: numerical - planets and satellites: formation - protoplanetary disks", '1. Introduction': "Planets form in gaseous disks around young stars. Initially, micron-sized dust grains are present in such a disk. It is generally thought that planets form from planetesimals of the order of a kilometer in size, which are large enough for gravity to play an important role (Safronov 1969; Pollack et al. 1996; Benz 2000). However, the first step in the planet formation process, the coagulation of dust grains into planetesimals, is still an unsolved problem. Although micron-sized grains can quite easily coagulate to mm-sized particles (Dominik & Tielens 1997; Ormel et al. 2007), typical relative velocities quickly become so large that particles fragment or bounce rather than stick (Blum & Wurm 2000; Zsom et al. 2010). The typical size at which coagulation is no longer feasible is referred to as the fragmentation or bouncing barrier. Observations at radio wavelengths indeed infer a large reservoir of these mm-to-cm size particles (pebbles) ( e.g. Testi et al. 2014).  \nEven if it were possible to coagulate uninterruptedly from micron-sizes to large bodies, the growth process would face another challenge: that of radial drift. As soon as particles reach sizes at which they aerodynamically start to decouple from the gas, they lose angular momentum and spiral towards the star (Whipple 1972; Weidenschilling 1977). Radial drift is a fast process and hence poses a challenge to planetesimal formation theories, as planetesimals should form before the solids are lost to the star.  \nA promising mechanism to quickly form planetesimals directly from pebbles is streaming instability (Youdin & Goodman 2005; Johansen et al. 2007; Johansen & Youdin 2007). Streaming instability leads to clumping of pebbles that subsequently become gravitationally unstable and form planetesimals. In order for streaming instabilities to occur, the metallicity of the disk has to be locally enhanced above the typical value of 1% (Johansen et al. 2009; Bai & Stone 2010; Dra¸˙zkowska & Dullemond 2014). However, if one considers a smooth disk ( i.e. , without any radial pressure bumps or special locations), solids are flushed to the star in an inside-out manner due to dust growth and radial drift ( e.g. , Birnstiel et al. (2010a, 2012); Krijt et al. (2016b); Sato et al. (2016)) 1 . Radial drift lowers the metallicity of the disk and - at first sight - works against creating conditions favorable to streaming instability.  \nThe viability of a pile-up of solids at the snowline has been a topic of interest for considerable time. Several mechanisms have  \nA special location that can a GLYPH<11> ect the distribution of solids in the disk is the snowline. The snowline is the radius from the central star interior to which all water is in the form of vapor. Recent observations of structures in protoplanetary disks (for example HL Tau; ALMA Partnership et al. 2015) have been proposed to be associated with icelines (not necessarily water), e.g. , Banzatti et al. (2015); Okuzumi et al. (2016); Nomura et al. (2016).  \nbeen proposed, such as the creation of a pressure maximum ( e.g. Brauer et al. 2008) or a deadzone ( e.g. Kretke & Lin 2007). Another pile-up mechanism is water di GLYPH<11> usion and condensation: water vapor in the inner disk can di GLYPH<11> use outward across the snowline and condense, thereby increasing the solids reservoir just outside the snowline. Works that have investigated the role of water di GLYPH<11> usion and condensation on the growth and distribution of solids include Stevenson & Lunine (1988); Cuzzi & Zahnle (2004); Ciesla & Cuzzi (2006); Ros & Johansen (2013); Estrada et al. (2016); Armitage et al. (2016) and, recently, Krijt et al. (2016a), which focuses on the atmospheric snowline. In this paper, we investigate the role of the radial snowline to enhance the solids-to-gas ratio.  \nThe works most relevant to our problem are Stevenson & Lunine (1988) and Ros & Johansen (2013). Stevenson & Lunine (1988) investigated the role of outward di GLYPH<11> usion of water vapor across the snowline. In their work, a 'cold-finger' e GLYPH<11> ect was considered: water vapor from the inner disk condenses on (preexisting) solids beyond the snowline. In this way, they found a steep increase in the solid density beyond the snowline and hypothesized that this could accelerate the formation of Jupiter. However, dynamical aspects of the disk, such as radial drift of solids and accretion of gas, were not taken into account; except for the vapor, their model is essentially static. More recently, Ros & Johansen (2013) also simulated pebble growth by condensation, adopting a Monte Carlo method. In their simulation, the purely icy particles grow by condensation, until the point that they are 'lost' by evaporation (modelled as an instantaneous process). Because their work does not feature silicate 'seeds' that are immune to evaporation, the remaining ice particles can only grow in size. Consequently, a steady-state is never reached, also not because the simulation domain does not connect to the rest of the disk: there is no inflow of small pebbles or removal of vapor-rich gas.  \nIn this work we consider a model that emphasizes the dynamic nature of the disk, characterized by the accretion rates of pebbles and gas. Our model of the iceline therefore includes inflow of pebbles (at the outer boundary) and outflow of vapor (at the inner boundary) - a design that results in a steady-state. Our goal is to investigate whether-or-not this more dynamic design will result in a solids enhancement that could trigger streaming instabilities. In particular, we will investigate how the solids-togas ratio will depend on the properties of the disk: the accretion rates (of solids and gas), the turbulence strength, and the aerodynamical properties of the drifting pebbles. We will also investigate which side of the snowline the solids-to-gas ratio is boosted most: interior or exterior.  \nWe use a characteristic particle method, assuming a single pebble size at each time and location in the disk (Ormel 2014; Sato et al. 2016). We numerically solve the partial differential equations that govern the transport and condensation / evaporation processes for solids and vapor. We consider two model designs for the composition of pebbles: 'single seed' (all silicates in one core) and 'many seeds' (many tiny silicates grains). Our model also accounts for feedback of solids on the gas (a.k.a. collective e GLYPH<11> ects) and feedback of water vapor on the scale height of the disk. A semi-analytic model is presented that approximates the numerical results, enabling us to carry out parameter searches.  \nWe start out with a description of our model in Sect. 2. In Sect. 4 we discuss numerical results for a fiducial set of model parameters in detail. We first describe the time-dependent numerical results in Sect. 4.1, after which we focus on the steadystate solution in Sect. 4.2. In Sect. 5 we summarize the semi-  \nanalytical approximate model (which is described in more detail in Appendix C) and present our results concerning streaming instability conditions. We discuss our results in Sect. 6 and present our main conclusions in Sect. 7. A list of frequently used symbols can be found in Table A.1.", '2.1. Disk model': 'Throughout this paper, we use the GLYPH<11> -accretion model (Shakura & Sunyaev 1973) to model disks. We fix the mass of the central star M ? = M GLYPH<12> . We fix the temperature profile T ( r ) to:  \nT = 150 GLYPH<18> r 3 au GLYPH<19> GLYPH<0> 1 = 2 K (1)  \nwhere r is the radial distance from the star. This temperature power-law profile corresponds to a passively-irradiated disk and a solar-mass star (Kenyon & Hartmann 1987; Armitage et al. 2016); for simplicity, we neglect viscous heating. The isothermal sound speed cs is related to the temperature as:  \ncs = s kBT GLYPH<22> (2)  \nwith kB the Boltzmann constant and GLYPH<22> the molecular weight of the gas, which is 2.34 times the proton mass for a typical solar metallicity gas. The disk scale height H gas is given by:  \nH gas = cs GLYPH<10> GLYPH<25> 0 : 033 GLYPH<18> r 1 au GLYPH<19> 5 = 4 au (3)  \nwhere GLYPH<10> is the Keplerian orbital frequency. The gas accretes to the central star at a rate ˙ M gas. Given GLYPH<11> and ˙ M gas, the steady-state gas surface density is (Lynden-Bell & Pringle 1974):  \nGLYPH<6> gas = ˙ M gas 2 GLYPH<25> r v gas = ˙ M gas 3 GLYPH<25>GLYPH<23> (4)  \nwhere v gas = 3 GLYPH<23>= 2 r is the radial speed of the accreted gas, with GLYPH<23> the turbulent viscosity, which is given by:  \nGLYPH<23> = GLYPH<11> csH gas (5)  \nwhere GLYPH<11> is the dimensionless turbulence parameter. As a result of our choices for the temperature and scale height profiles, GLYPH<6> gas( r ) is a power-law with index GLYPH<0> 1. We take the turbulent gas di GLYPH<11> usivity D gas equal to the turbulent viscosity GLYPH<23> .', '2.2. Radial motion of ice, silicates and vapor': "In this study, we are interested in the pebble-sized fraction of the solids - we do not consider the smaller-sized dust. Pebbles are characterized by the fact that they are aerodynamically partly decoupled from the gas disk. The gas disk is partly pressuresupported, and therefore rotates at a sub-Keplerian velocity. The pebbles do not feel the pressure gradient, and therefore tend to move toward Keplerian orbital velocities. This means that pebbles are moving through a more slowly rotating gas disk, thereby losing angular momentum through gas drag. The loss of angular momentum results in pebbles spiralling ('drifting') towards the central star. The extent to which pebbles drift depends on the stopping time t stop, which measures how strongly the pebbles are coupled to the gas - it is the time after which any initial momentum is lost due to gas friction.  \nConcerning the stopping times of pebbles we can distinguish between two regimes: in the Stokes regime, the particle size is larger than the mean-free path of the gas molecules l mfp, and the stopping time is calculated in a fluid description; in the Epstein regime, the particle size is smaller than the mean-free path of the gas molecules, and a particle description is needed instead. The mean-free path l mfp is given by:  \nl mfp = GLYPH<22> p 2 GLYPH<26> gas GLYPH<27> mol (6)  \nwhere GLYPH<27> mol is the molecular collision cross-section, and GLYPH<26> gas is the gas density. We take GLYPH<27> mol = 2 GLYPH<2> 10 GLYPH<0> 15 cm 2 as the collisional cross-section of molecular hydrogen (Chapman & Cowling 1970). The stopping time is given by:  \nt stop = 8 > > > > > > > < > > > > > > > : GLYPH<26> GLYPH<15> ; p s p v th GLYPH<26> gas ; (Epstein: s p < 9 4 l mfp) 4 GLYPH<26> GLYPH<15> ; p s 2 p 9 v th GLYPH<26> gas ; (Stokes: s p > 9 4 l mfp) (7)  \nwhere s p is the particle radius, GLYPH<26> GLYPH<15> ; p is the particle internal density, and v th is the thermal velocity of the gas molecules, defined as v th = p 8 =GLYPH<25> cs . The dimensionless stopping time (sometimes referred to as the Stokes number) is GLYPH<28> S = t stop GLYPH<10> .  \nThe radial drift velocity of pebbles v peb is given by (Weidenschilling 1977; Nakagawa et al. 1986):  \nv peb = v gas + 2 GLYPH<17>v K GLYPH<28> S 1 + GLYPH<28> 2 S (8)  \nwhere GLYPH<17>v K is the magnitude of the azimuthal motion of the gas disk below the Keplerian velocity v K :  \nGLYPH<17>v K = GLYPH<0> 1 2 c 2 s v K @ log P @ log r (9)  \nwith P the gas pressure. Under normal disk conditions, the solids-to-gas ratio is about 1 in 100, and it is reasonable to neglect the back-reaction of the solids on the gas. In this paper we are interested in situations where the solids-to-gas ratio is boosted, however, and therefore we study the e GLYPH<11> ects of the backreaction of solids on the gas, or 'collective e GLYPH<11> ects', as well. Collective e GLYPH<11> ects have been calculated by Nakagawa et al. (1986). However, their formula was derived in the inviscid limit. In Appendix B we have calculated the drift velocity, accounting for both the backreaction and viscous forces. When including the back-reaction, we replace Eq. (8) by Eq. (B.7).  \nWater vapor (denoted by the subscript Z ) mixes with the gas and gets accreted to the central star at the same velocity as the gas.", '2.3. Evaporation and condensation': 'The rate of change in mass ( dm p = dt ) of a spherical icy pebble is given by (Lichtenegger & Komle 1991; Ciesla & Cuzzi 2006; Ros & Johansen 2013):  \ndm p dt = 4 GLYPH<25> s 2 p v th ; Z GLYPH<26> Z 1 GLYPH<0> P eq PZ ! (10)  \nwhere s p is the radius of the pebble, v th ; Z is the thermal velocity of vapor particles, GLYPH<26> Z is the vapor density, PZ is the vapor pressure, and P eq is the saturated or equilibrium pressure, which is given by the Clausius-Clapeyron equation:  \nP eq = P eq ; 0 e GLYPH<0> Ta = T (11)  \nwhere Ta and P eq ; 0 are constants depending on the species. For water, Ta = 6062 K and P eq ; 0 = 1 : 14 GLYPH<2> 10 13 g cm GLYPH<0> 1 s GLYPH<0> 2 (Lichtenegger & Komle 1991).  \nAssuming spherical pebbles and one typical particle mass m p at each location and time in the disk ( i.e. , the particle size distribution is narrow), and using the ideal gas law to write PZ in terms of the (midplane) vapor density GLYPH<26> Z , we can rewrite Eq. (10) into a change in the ice surface density GLYPH<6> ice due to condensation and evaporation:  \n˙ GLYPH<6> ice ; C = E = R c GLYPH<6> ice GLYPH<6> Z GLYPH<0> R e GLYPH<6> ice (12)  \nwhere the dot denotes the time derivative. In Eq. (12) the vapor surface density GLYPH<6> Z is related to the midplane vapor density as GLYPH<26> Z = GLYPH<6> Z = p 2 GLYPH<25> H gas and R c and R e are the condensation rate and evaporation rate, respectively, defined as:  \nRc = 8 s kBT GLYPH<22> Z s 2 p m p H gas (13)  \nRe = 8 p 2 GLYPH<25> s 2 p m p r GLYPH<22> Z kBT P eq (14)  \nwith kB the Boltzmann constant, and GLYPH<22> Z = 18 m H the mean molecular weight of water vapor, where m H is the proton mass. Evaporating ice turns into vapor, and condensing vapor turns into ice, and therefore we find for the vapor source terms:  \n˙ GLYPH<6> Z ; C = E = GLYPH<0> ˙ GLYPH<6> ice ; C = E = GLYPH<0> R c GLYPH<6> ice GLYPH<6> Z + R e GLYPH<6> ice (15)', '2.4. Internal structure and composition of pebbles': "Weassume that in the outer disk, half of the mass in pebbles consists of water ice, and the other half consists of silicates (Lodders 2003; Morbidelli et al. 2015). That is, the dust-fraction in pebbles in the outer disk GLYPH<16> = 0.5. Adopting GLYPH<26> GLYPH<15> ; ice = 1 g cm GLYPH<0> 3 as the density of a pure ice particle and GLYPH<26> GLYPH<15> ; sil = 3 g cm GLYPH<0> 3 as the density of a pure silicate particle, we find that in the outer disk, the internal density of a pebble is GLYPH<26> GLYPH<15> ; p = 1 : 5 g cm GLYPH<0> 3 . Generally, the internal density of an ice / silicate pebble is given by:  \nGLYPH<26> GLYPH<15> ; p = GLYPH<26> GLYPH<15> ; ice GLYPH<26> GLYPH<15> ; sil( m sil + m ice) m sil GLYPH<26> GLYPH<15> ; ice + m ice GLYPH<26> GLYPH<15> ; sil (16)  \nWe consider two possibilities for the internal structure of icy pebbles (Fig. 1).  \nwhere m ice is the mass in ice within the pebble and m sil is the mass in silicate within the pebble. Equation (10) is valid for a pure ice particle, while we use it for icy pebbles with (a) silicate core(s). We checked that using m ice instead of m p in the condensation and evaporation rates does not change the results significantly.  \n- - Single-seed model . In this model we view pebbles in the outer disk as balls of ice with a single core made of silicates. When the ice in such a pebble evaporates, the silicate core is left behind. In this framework, the number of pebbles does not change when crossing the snowline.\n- - Many-seeds model . In this model we view icy pebbles as lots of small silicate particles that are 'glued' together by ice. The small silicate particles are released by evaporation, and free silicates can stick onto icy pebbles. We do not 'resolve' the internal structure of the pebbles, and only trace the total amount of silicates 'locked' into icy pebbles. Following Saito & Sirono (2011), we assume that the released silicate particles are of micron size. Because they are so small, the silicates do not undergo any aerodynamic drift and follow the radial motion of the gas.  \nFig. 1. Schematic showing the di GLYPH<11> erence between the single-seed model and the many-seeds model for the internal structure of icy pebbles in the outer disk. In the single-seed model, the water ice layer (white) on the silicate core (red) becomes smaller at the evaporation front, and after complete evaporation of the ice, only one single silicate core remains. Since we take GLYPH<16> = 0 : 5, the silicate core contains half the mass of the original pebble: m core = m ice = 0 : 5 m p ; start. In the many-seeds model, pebbles in the outer disk consist of many micron-sized silicate grains (red) 'glued' together by water ice (white). As the ice evaporates, silicate grains are released as well. After complete evaporation of the ice, many micron-sized silicate grains remain. Again, with GLYPH<16> = 0 : 5, the total mass in silicate grains in a pebble in the outer disk is half of the total pebble mass, but this mass is now divided among many silicate grains that each have mass m sil. See text for more details.  \n<!-- image -->", '2.5. Transport equations': 'In this subsection we provide the systems of transport equations corresponding to both pebble composition models presented in Sect. 2.4. Our model is 1D; we integrate over the vertical dimension of the disk and follow column (surface) densities. Midplane densities are calculated assuming a Gaussian distribution in height with scaleheights H gas (for the gas and the vapor) and H peb (for the pebbles), respectively.', '2.5.1. Single-seed pebble model': 'In the single-seed model, we follow the vapor surface density GLYPH<6> Z , the ice surface density GLYPH<6> ice, and the number density of particles N p. The equations governing the time evolution of these quantities are:  \n@ GLYPH<6> ice @ t + r GLYPH<1> M ice = ˙ GLYPH<6> ice ; C = E (17)  \n@ GLYPH<6> Z @ t + r GLYPH<1> M Z = ˙ GLYPH<6> Z ; C = E (18)  \n@ N p @ t + r GLYPH<1> N p = 0 (19)  \nwhere the source terms ˙ GLYPH<6> ice ; C = E and ˙ GLYPH<6> Z ; C = E are defined in Eq. (12) and Eq. (15), M ice and M Z are the ice and vapor mass flux, respectively, and N p is the number flux of solid particles. These  \nfluxes are given by:  \nM ice = GLYPH<0> GLYPH<6> ice v peb GLYPH<0> D p GLYPH<6> gas r GLYPH<6> ice GLYPH<6> gas (20)  \nM Z = GLYPH<0> GLYPH<6> Z v gas GLYPH<0> D gas GLYPH<6> gas r GLYPH<6> Z GLYPH<6> gas (21)  \nN p = GLYPH<0> N p v peb GLYPH<0> D p N gas r N p N gas (22)  \nwhere v peb is the drift velocity of the pebbles given by Eq. (8) or by Eq. (B.7) when including collective e GLYPH<11> ects, N gas is the number of gas particles, and D p is the particle di GLYPH<11> usivity, which is related to D gas through the Schmidt number Sc (Youdin & Lithwick 2007):  \nD p = D gas GLYPH<2> Sc = Dgas 1 + GLYPH<28> 2 S (23)  \nSince in our models the stopping times are much smaller than unity, we set D p equal to D gas.  \nNote that a negative flux means that material is transported inwards, towards the star. Under the single-size approximation, the typical pebble mass m p is given by:  \nm p = m core + GLYPH<6> ice N p (24)  \nwith m core the mass of the bare silicate cores of the pebbles.', '2.5.2. Many-seeds pebble model': 'We implement the many-seeds model as follows. We now consider two additional populations of solids: dirt and free silicates. Dirt refers to the silicates that are captured in icy pebbles, whereas free silicates are small silicate grains without any ice coating. The transport equations for GLYPH<6> ice, GLYPH<6> Z , and N p are the same as in the single-seed pebble model (Eq. (17)-(19)). The transport equations for the dirt surface density GLYPH<6> dirt and free silicates surface density GLYPH<6> sil are given by:  \n@ GLYPH<6> dirt @ t + r GLYPH<1> M dirt = GLYPH<0> max[0 ; ( R e GLYPH<0> R c GLYPH<6> Z )] GLYPH<6> dirt + R s GLYPH<6> sil (25)  \n@ GLYPH<6> sil @ t + r GLYPH<1> M sil = max[0 ; ( R e GLYPH<0> R c GLYPH<6> Z )] GLYPH<6> dirt GLYPH<0> R s GLYPH<6> sil (26)  \nwhere M dirt and M sil are the dirt and free silicate mass fluxes, respectively. Free silicates are released from the pebbles at a rate proportional to the decrease of ice surface density due to evaporation and to the amount of dirt that can evaporate from within the pebbles: ( Re GLYPH<0> Rc GLYPH<6> vap) GLYPH<6> dirt, but only when there is net evaporation of ice. Silicate grains stick to icy pebbles at a collision / sticking rate R s:  \nR s = GLYPH<6> ice + GLYPH<6> dirt p 2 GLYPH<25> H sil m p GLYPH<1> GLYPH<1> v sil ; peb GLYPH<1> GLYPH<25> s 2 p (27)  \nwhere GLYPH<1> v sil ; peb is the relative velocity between the pebbles and the silicates and H sil is the scale height of the free silicates. Since the inward drift velocity of the pebbles is much larger than the  \noutward di GLYPH<11> usion velocity of the silicates, we take GLYPH<1> v sil ; peb GLYPH<25> v peb 2 .  \nThe dirt and silicate mass fluxes, M dirt and M sil, are given by:  \nM dirt = GLYPH<0> GLYPH<6> dirt v peb GLYPH<0> D p GLYPH<6> gas r GLYPH<6> dirt GLYPH<6> gas (28)  \nM sil = GLYPH<0> GLYPH<6> sil v gas GLYPH<0> D gas GLYPH<6> gas r GLYPH<6> sil GLYPH<6> gas (29)  \nwhere v peb is again the radial drift velocity of the pebbles, and silicates are advected with the gas.', '2.6. Effects of variable GLYPH<22>': 'When icy pebbles evaporate, the gas becomes enriched in water vapor. This means that the mean molecular weight GLYPH<22> of the gas increases:  \n1 GLYPH<22> = f H2O GLYPH<22> H2O + f gas GLYPH<22> gas (30)  \nwhere f H2O is the mass fraction of water molecules, GLYPH<22> H2O = 18 m H is the mean molecular weight of water, and f gas is the mass fraction of gas molecules (without the vapor contribution), with GLYPH<22> gas = 2 : 34 m H.  \nSince the sound speed cs is proportional to p 1 =GLYPH<22> , the sound speed decreases as more water vapor is added. The disk scale height H gas is proportional to cs and will therefore decrease as well: the disk becomes more vertically compact with increasing water vapor. An increase in GLYPH<22> across the water evaporation front has an e GLYPH<11> ect on the pebble drift velocity v peb (Eq. (8)) in this region, since v peb is proportional to 1 =GLYPH<22> through its dependence on cs , and sensitive to variations in the gas pressure P through @ log P =@ log r . The total gas pressure (hydrogen / helium plus water) is given by:  \nP = P gas + PZ = GLYPH<6> gas kBT p 2 GLYPH<25> H gas GLYPH<22> H2 = He + GLYPH<6> ZkBT p 2 GLYPH<25> H gas GLYPH<22> H2O (31)  \nWe take the viscosity GLYPH<23> to be independent of variations in GLYPH<22> . Therefore, GLYPH<6> gas = ˙ M gas = 3 GLYPH<25>GLYPH<23> is independent of GLYPH<22> as well. However, the H2O-enriched gas is characterized by a reduced scaleheight, since H gas / p 1 =GLYPH<22> , increasing the midplane gas pressure P gas as P gas / p GLYPH<22> . Evaporation hence increases the pressure gradient and the headwind of the gas through Eq. (9): pebbles tend to drift faster. Finally, the drift velocity v peb also depends on the stopping time GLYPH<28> S , which in turn depends on the gas density and the thermal velocity - both quantities that depend on GLYPH<22> . We take these dependencies into account as well.  \nWe now define two models that we study in the following sections. In the simple model, we do not include collective effects and the variability of the mean molecular weight GLYPH<22> . In the simple model, the sound speed and scale height are independent of variations in GLYPH<22> , GLYPH<17>v K (Eq. (9)) is constant and the midplane gas density is GLYPH<22> -independent: GLYPH<26> gas = GLYPH<6> gas = p 2 GLYPH<25> H gas. In the complete model, we do take collective e GLYPH<11> ects and variations of GLYPH<22> into account, as described above.', '2.7. Input parameters': "The four key input parameters that can be varied in our model are:  \n- -The gas accretion rate ˙ M gas.\n- -The ratio between the pebble accretion rate and the gas accretion rate, denoted by the symbol F s =g , and defined as:  \nF s =g GLYPH<17> ˙ M peb = ˙ M gas : (32)  \n- -The dimensionless turbulence parameter GLYPH<11> .\n- -The initial stopping time of pebbles just outside the snowline (at 3 au) denoted by the symbol GLYPH<28> 3. This parameter is simply a proxy for the mass of the incoming pebbles.  \nTogether, these input parameters define the 'zero-model' of the disk: the gas and solids distribution in the case of no condensation and evaporation. The gas surface density profile GLYPH<6> gas is defined through Eq. (4), where the turbulent viscosity GLYPH<23> is determined by GLYPH<11> through Eq. (5). With GLYPH<6> gas and the stopping time normalization GLYPH<28> 3 at 3 au, we can find GLYPH<28> ( r ) at any other radius r . We can then also determine the pebbles surface density GLYPH<6> peb for the zero-model: it is given by F s =g ˙ M gas = 2 GLYPH<25> r v peb( r ), where the pebble velocity v peb( r ) depends on GLYPH<28> ( r ) and is given by Eq. (8) (without backreaction) or Eq. (B.7) (including backreaction).", '2.8. Model assumptions': "In constructing our model we employed several assumptions, which we discuss below.  \n- -We do not include a particle size distribution, but rather consider one typical pebble size at each time and location in the disk. In the many-seeds model, we do distinguish between the population of small silicate grains and large icy pebbles, however. We come back to the importance of condensation onto small grains rather than onto pebbles in the discussion (Sect. 6).\n- -Our model is 1D, and has a straight vertical snowline. In reality, the vertical snowline is curved rather than straight, because the temperature varies with height above the midplane ( e.g. Oka et al. (2011)). However, Ros & Johansen (2013) have shown that particle growth due to di GLYPH<11> usion and subsequent condensation of water across the vertical snowline is negligible compared to the same process across the radial snowline.\n- -Our model does not account for coagulation among pebbles: the enhancement of the solids-to-gas ratio is solely due to di GLYPH<11> usion and condensation. We remain agnostic as to what happens to the pebbles in terms of coagulation and fragmentation until they reach the evaporation front. In the manyseeds model, we do not take into account coagulation of silicate particles. We will come back to this point in the discussion (Sect. 6).\n- -Similarly, we assume that the scale height of the vapor - and, in the many-seeds implementation, of the free silicates - is always equal to the scale height of the disk; in other words, when an icy pebble evaporates, the released water vapor immediately distributes itself over the same vertical extent as the background gas, whereas it originates from the icy pebbles that reside in a thinner disk, characterized by a scale height H peb:  \nH peb = H gas r GLYPH<11> GLYPH<28> S + GLYPH<11> : (33)  \nThe assumption of rapid vertical mixing extends to the calculation of the mean molecular weight of the gas in the complete model. Also there, it is assumed that the water vapor is distributed evenly throughout the gas at all times. In Sect. 4.2, we verify that the assumption of taking the vapor scale height equal to the disk scale height rather than the pebbles scale height, does not have a large e GLYPH<11> ect on the results.  \n- -Weassume a steady-state gas distribution, that is, ˙ M gas and T remain constant. This assumption is justified since the solidsto-gas enhancement takes place on a radial scale GLYPH<1> r that is small compared to the size of the disk. In other words, the local viscous evolution timescale ( GLYPH<1> r 2 =GLYPH<23> ), on which the effect happens, is much smaller than the global viscous evolution timescale of the disk, on which the gas accretion rate and global temperature profile change. The temperature profile (Eq. (1)) corresponds to a passively-irradiated disk. However, in case of viscous heating T ( r ) will still be described by a power-law (Armitage 2015). We have checked the e GLYPH<11> ect of adopting a viscous heating temperature profile ( T ( r ) / r GLYPH<0> 3 = 4 ) in our numerical simulation. The key consequence is that the location of the snowline is shifted in radius; all other model results are very similar.\n- -We adopt an GLYPH<11> -viscosity model, with constant GLYPH<11> (the value of which we allow to vary between very high and very low values), independent of the solids-to-gas ratio. However, if turbulence is driven by the magneto-rotational instability, the local turbulence strength at the midplane depends on the local grain size and abundance. So-called 'dead zones', where turbulence is strongly reduced, could occur at locations in the disk with a sharp increase in the solids-to-gas ratio (Kretke & Lin 2007). Since our results depend strongly on the amount of turbulence near the snowline, the creation of dead zones would change our results.\n- -Similarly, we do not consider the interplay between (small) grains and the local temperature: we do not solve for the temperature profile self-consistently, but take it as constant, given by Eq. (1). However, in the many-seeds model, a large amount of small silicate grains is released in the evaporation front. This means that locally the opacity could increase. Conceivably, this could lead to a steeper temperature profile in the region near the snowline and therefore to a more narrow evaporation front. Including these e GLYPH<11> ects requires a radiation transfer model, which is beyond the scope of this work.\n- -In the many-seeds implementation, we do not include condensation onto the small, free silicate particles. We also do not model the porosity of pebbles. We come back to these issues in Sect. 6.3.", '3. Numerical implementation': 'We make use of the open source partial di GLYPH<11> erential equations solver FiPy (Guyer et al. 2009) to solve the transport equations numerically. The coupled systems of equations are solved on a cylindrical 1D grid that ranges from 1 to 5 au.  \nWe implement boundary conditions as follows. At the outer boundary r 2 = 5 au, the pebble mass flux F s = g ˙ M gas is fixed 3 . We convert the input parameter GLYPH<28> 3 to a physical pebble size s p ; start using Eq. (7), and since our model does not include coagulation,  \nFig. 3. Rate of change of the ice surface density at the peak location, as function of time. The rate of change increases until the icy pebble front has reached the peak location, after which the peak continues to grow at an exponentially decreasing rate, while the system is converging to its steady-state solution.  \n<!-- image -->  \nwe set the pebble size (which can be converted to a pebble mass m p ; start) at the outer boundary r 2 to this value. For the singleseed model, the mass of the bare silicate core is then given by m core = GLYPH<16> m p ; start = m p ; start = 2. At the inner boundary r 1 = 1 au, we implement outflow boundary conditions: solid particles and water vapor are accreted to the central star at a constant rate. We also demand that GLYPH<6> ice = 0 at r 1 and GLYPH<6> Z = 0 at r 2, for obvious reasons. Additionally, in the many-seeds implementation we have GLYPH<6> dirt = GLYPH<6> ice = 0 at r 1 and GLYPH<6> sil = GLYPH<6> Z = 0 at r 2.', '3.1. Time-dependent': 'In our time-dependent method, we adopt small time steps ( GLYPH<1> t GLYPH<24> 100 yr) and a small value for the allowed residuals, to ensure mass conservation by demanding the residuals 4 to be small. This method requires a lot of computational power because the system of partial di GLYPH<11> erential equations is numerically very sti GLYPH<11> , but the advantage is that we find not only the steady-state solution to the equations, but are able to follow the solution in time. This allows us to get a measure of the time needed for the solution to settle into a steady-state. We present the time-dependent results for the simple, single-seed, model with fiducial input parameters in Sect. 4.1.', '3.2. Time-independent': "Our time-independent method solves for the steady-state solution directly, by taking very large time-steps (in total 40 time steps that increase to 10 7 yr) and decreasing the residuals threshold after each run, taking the solution from the previous run as the new input. We continue this process until we reach convergence. We present the steady-state results for the 'simple, single-seed', 'simple, many-seeds', 'complete, single-seed', and 'complete, many-seeds' models with fiducial input parameters in Sect. 4.2.  \nFig. 2. The solution to the transport equations at di GLYPH<11> erent points in time, for the fiducial model parameters listed in Table 1. At t = 10 5 yr, approximately 90% of the steady-state peak in the ice surface density has formed. (a) Surface densities GLYPH<6> of ice (solid lines) and vapor (dashed lines). The dotted line corresponds to the steady-state advection-only vapor surface density profile. (b) Midplane pebbles-to-gas ratio GLYPH<26> peb =GLYPH<26> gas. (c) Typical pebble mass m p. (d) Typical pebble internal density GLYPH<26> GLYPH<15> ; p.  \n<!-- image -->  \nTable 1. Fiducial model parameters  \nGLYPH<2>  \nNotes. See Sect. 2.7. GLYPH<28> 3 translates into a physical pebble size at 3 au for the zero-model.", '4. Results - Fiducial model': 'In this section we introduce a standard set of parameters which we then use to illustrate di GLYPH<11> erent aspects of the model. The fiducial model is defined by the parameters listed in Table 1. We first discuss the time-dependent results for the simple, singleseed, fiducial model in Sect. 4.1. We then go on to discuss the steady-state results for the other variants of the fiducial model in Sect. 4.2.', '4.1. Time-dependent solution': "In Fig. 2 we plot several quantities as function of radial distance from the star r , at four points in time. In Fig. 2a, the surface densities of ice (solid lines) and vapor (dashed lines) are plotted. Initially ( t = 0), the ice and vapor surface density are zero, GLYPH<6> ice = GLYPH<6> Z = 0. Icy pebbles are drifting in from the outer disk (right) to the inner disk (left). At the point where ˙ GLYPH<6> Z ; C = E > 0 the icy component of the pebbles evaporate. At the snowline r snow GLYPH<25> 2 : 1 au all the H2O is in the gas phase. Since the vapor is advected with the gas velocity v gas, which is much smaller than the pebble velocity v peb, the vapor surface density (dashed line) quickly exceeds that of the ice (solid line). The vapor density increases until the steady-state vapor distribution GLYPH<6> Z ; a has been reached. GLYPH<6> Z ; a is given by 5 :  \nGLYPH<6> Z ; a = ˙ M gas 3 GLYPH<25>GLYPH<23> (34)  \nand is plotted by the dotted curve. Meanwhile, due to the outward di GLYPH<11> usion of water vapor across the snowline and subse-  \nFig. 4. Steady-state results for the fiducial parameters, for the simple model with single-seed evaporation (upper two panels) and many-seeds evaporation (lower two panels). (a) Surface densities of ice, silicates and vapor ( GLYPH<6> ice, GLYPH<6> sil, and GLYPH<6> Z , respectively) for the single-seed model. The silicates are locked in the icy pebbles outside the snowline, whereas they are bare grains interior to the snowline. (b) Midplane pebbles-to-gas ratio GLYPH<26> peb =GLYPH<26> gas, and pebbles-to-gas column density ratio GLYPH<6> peb = GLYPH<6> gas for the single-seed model. (c) Same as panel (a) but for the many-seeds model. The silicates are now divided in two populations: silicates that are locked up in icy pebbles ( GLYPH<6> dirt), and free silicates ( GLYPH<6> sil). Note that in panels (a) and (c) the x -axis ranges between r = 1 : 8 GLYPH<0> 2 : 8 au, for clarity. (d) Same as (b), but for the many-seeds model and with an additional line for the midplane 'free silicates'-to-gas ratio GLYPH<26> sil =GLYPH<26> gas.  \n<!-- image -->  \nquent condensation onto the inward-drifting pebbles, the ice surface density just exterior to the snowline increases. Di GLYPH<11> usion of vapor therefore results in a distinct bump in the ice profile and we denote the location corresponding to maximum GLYPH<6> ice the peak radius r peak. A similar maximum is seen in the solids-to-gas ratio (Fig. 2b), which rises to GLYPH<24> 0 : 45. After t = 2 GLYPH<2> 10 5 yr a steady state has been reached.  \nThe enhancement of the pebbles-to-gas ratio is reflected by the typical pebble mass m p and density GLYPH<26> GLYPH<15> ; p, plotted in Fig. 2c and Fig. 2d, respectively. Pebbles at the inner boundary are half as massive as pebbles at the outer boundary, because at the inner boundary all the ice has evaporated o GLYPH<11> the pebbles, and the dust fraction of pebbles GLYPH<16> = 0 : 5. The typical pebble mass just outside the snowline increases over time, indicating that water vapor condenses onto the pebbles. This is also shown by the behaviour of the typical pebble density GLYPH<26> GLYPH<15> ; p. At the outer boundary, GLYPH<26> GLYPH<15> ; p = 1 : 5 g cm GLYPH<0> 3 , corresponding to pebbles containing as much water ice as silicates in mass; at the inner boundary, GLYPH<26> GLYPH<15> ; p = 3 g cm GLYPH<0> 3 , corresponding to pure silicate pebbles; whereas  \njust outside the snowline, GLYPH<26> GLYPH<15> ; p decreases over time to values close to unity - corresponding to pure water ice pebbles. A decrease in particle internal density outside the snowline due to water condensation was also found by Estrada et al. (2016).  \nIn Fig. 3 we plot the growth rate of the peak in GLYPH<6> ice outside the snowline as function of time t . At the start of the simulation, ice is not yet present at the location of the peak. After the icy pebbles have reached the location of the eventual peak, at t GLYPH<24> 2 GLYPH<2> 10 4 yr, the rate of growth decreases exponentially - the peak continues to grow, but at an ever slower rate, while the system is converging to its steady-state solution. The e -folding growth timescale is GLYPH<24> 3 GLYPH<2> 10 4 yr, which is approximately equal to the gas advection timescale ( r peak GLYPH<0> r snow) =v gas; the time it takes for the vapor to traverse the peak and contribute to the steep vapor concentration gradient along which outward di GLYPH<11> usion can take place. We will come back to the peak formation timescale in Sect. 5. The time in which 90% of the peak forms is GLYPH<24> 10 5 yr. Note that we have assumed that the incoming pebble mass flux does not change over time ( F s =g is constant at the outer boundary r 2). If the pebble flux would decrease on a timescale shorter than the  \ntime it takes to form the peak, a steady-state would not exist. If the pebble flux would decrease on a timescale longer than it takes to form the peak, we would reach a quasi-steady-state solution; in that case our model is fully applicable.", '4.2. Steady-state solution': 'The time-dependent solution described in Sect. 4.1 converges to the steady-state solution. In this section we take a closer look at the steady-state solution for the fiducial parameters, for the model designs discussed in Sect. 2.', '4.2.1. Single-seed versus many-seeds': 'In Fig. 4 we compare the simple single-seed model (upper two panels) with the simple many-seeds model (lower two panels). The left panels show the surface densities of ice, vapor, and in the many-seeds case, also of locked silicates (dirt) and free silicates. The ice profiles are very similar in both models, but in the many-seeds case there is an additional peak in the dirt surface density. This is because the free silicates behave like vapor: they are released throughout the evaporation front, and follow the gas profile interior to the snowline. Turbulent di GLYPH<11> usion not only causes vapor to be transported outward and condense onto the pebbles, but leads to the same e GLYPH<11> ect for the small silicates - they di GLYPH<11> use outward and stick onto the pebbles, resulting in a peak in the dirt surface density distribution.  \nIn the right panels we plot pebbles-to-gas midplane density ratios GLYPH<26> peb =GLYPH<26> gas and column density ratios GLYPH<6> peb = GLYPH<6> gas. In the manyseeds model, the pebbles-to-gas ratio is the sum of the dirt-togas ratio and the ice-to-gas ratio. In both models, a clear peak is present just outside the snowline. The midplane density ratio is larger than the column density ratio due to settling of pebbles. In the many-seeds case, we also plot the free silicates-to-gas midplane density ratio GLYPH<26> sil =GLYPH<26> gas. This ratio is constant interior to the snowline because the silicates are perfectly coupled to the gas. In the many-seeds model, small silicate grains di GLYPH<11> use outward across the snowline, whereas in the single-seed model, there is no significant outward di GLYPH<11> usion of the larger silicate pebbles due to the absence of a steep concentration gradient. Therefore, the peak in the pebbles-to-gas ratio outside the snowline is larger in the many-seeds model than in the single-seed model.', '4.2.2. Simple versus complete models': 'Figure 5 shows the results for the complete model with the single-seed pebble design (upper two panels, a and b) and with the many-seeds pebble design (lower two panels, c and d). In Fig. 5a (single-seed) and Fig. 5c (many-seeds), we plot the surface densities of ice, vapor, and in the many-seeds case the surface densities of locked silicates (dirt) and free silicates (sil). Dashed lines denote the simple model results of Fig. 4a and Fig. 4c, for comparison. The grey lines correspond to the mean molecular weight of the gas. From the outer disk to the inner disk, the mean molecular weight increases from 2.34 m H to 3.11 m H across the snowline. Clearly, the peak in the ice surface density is higher and broader for the complete model than for the simple model. This is because collective e GLYPH<11> ects (which reduce the pebble drift velocity) outweigh the e GLYPH<11> ect of the enhanced gas pressure gradient in the evaporation front (which enhances the pebble drift velocity). In the complete model, therefore, pebbles e GLYPH<11> ectively have a smaller radial velocity than in the simple model, and therefore the resulting ice peak is higher and  \nbroader. This is also the case in the many-seeds implementation. The width of the peak strengthens the justification of the vertical mixing assumption (see below).  \nIn Fig. 5b and Fig. 5d we compare the midplane pebbles-togas ratios between the simple and the complete model. Collective e GLYPH<11> ects help to boost the pebbles-to-gas ratio, for the reason outlined above.', '4.2.3. Validation of vertical mixing assumption': "As discussed in Sect. 2.8, we assume rapid vertical mixing of vapor. However, it could be argued that (part of) the released vapor recondenses onto icy pebbles before relaxing to the same vertical distribution as the gas, thereby boosting the condensation rate. This assumption also extends to the small silicate grains in the many-seeds model. Similarly, one could argue that the silicate grains can collide with a pebble before they can di GLYPH<11> use to higher vertical layers ( e.g. Krijt & Ciesla (2016)). In order to investigate what the di GLYPH<11> erence is in terms of the steady-state ice distribution between two extreme cases, we take the simple, single-seed, fiducial model and run it once with HZ = H gas and once with HZ = H peb. We compare the steady-state outcomes of the two cases in Fig. 6. Even though the location of the snowline and of the pebbles-to-gas peak di GLYPH<11> ers a bit between the two cases, the height of the peak in the midplane pebbles-to-gas ratio increases only by GLYPH<24> 0 : 3% if one takes HZ = H peb instead of HZ = H gas. Similarly, assuming a smaller scale height for the silicate grains in the evaporation region would probably not significantly change the height of the peak in the 'locked' silicates, since the silicate grains behave like vapor. In the rest of this paper, we assume HZ = H sil = H gas.", '5. Streaming instability conditions': 'Carrera et al. (2015) have investigated for what values of the solids-to-gas column density ratio streaming instability occurs, as function of the stopping time of the solid particles. They found that, typically, the solids-to-gas column density ratio should exceed GLYPH<24> 2% for stopping times GLYPH<24> 5 GLYPH<2> 10 GLYPH<0> 2 . However, these results should be interpreted with caution, since Carrera et al. (2015) only accounted for self-driven (Kelvin-Helmholtz) turbulence. We do include global turbulence, and therefore cannot simply use their conditions for streaming instability. A more robust condition for streaming instability in the presence of global turbulence is probably a midplane solids-to-gas ratio that reaches order unity (Johansen & Youdin 2007).  \nIn this section we investigate under what disk conditions, the midplane solids-to-gas ratio near the snowline gets most enhanced due to the e GLYPH<11> ect of water di GLYPH<11> usion and condensation. To this end, we have constructed a semi-analytic approximate model to be able to quickly test for many di GLYPH<11> erent values of the input parameters (Sect. 2.7), within the single-seed implementation (Sect. 2.4). Before presenting the results, we will provide a short summary of the semi-analytical model, which is discussed in greater detail in Appendix C.', '5.1. Semi-analytical model: Summary': "Our semi-analytical model is an analytic prescription to find the approximate location, height and width of the ice surface density at the location of the ice peak, r peak, of the steady-state (timeindependent) solution (see Sect. 3.2 for our numerical method to find the steady-state solution to the transport equations). For sim-  \nFig. 5. Comparison of steady-state results for the fiducial parameters between the simple model (dashed lines) and the complete model (solid lines), which includes collective e GLYPH<11> ects and the e GLYPH<11> ects of the variation of the mean molecular weight GLYPH<22> . The upper two panels correspond to the single-seed evaporation model; the lower two panels correspond to the many-seeds evaporation model. In the latter case, we made the conservative choice of not including the free silicates in the back-reaction onto the gas. (a) Surface densities of ice ( GLYPH<6> ice), vapor ( GLYPH<6> Z ) and silicates ( GLYPH<6> sil) for the single-seed model. The grey line indicates the mean molecular weight of the gas in units of proton mass m H. (b) Midplane pebbles-to-gas ratio GLYPH<26> peb =GLYPH<26> gas. (c) Same as (a), but the silicates are now divided into 'free silicates' ( GLYPH<6> sil) and 'locked silicates' ( GLYPH<6> dirt). For clarity, we only show the clean GLYPH<6> ice profile for comparison. (d) Same as (b), with an additional line denoting the 'free silicates'-to-gas ratio GLYPH<26> sil =GLYPH<26> gas.  \n<!-- image -->  \nplicity, we have tailored the semi-analytical model towards the 'single-seed' design. (An extension to also include the 'manyseeds' design will be considered in a future work).  \nThe model works by iteratively solving a set of analytical equations, until convergence is reached. The key parameter is the ratio between the pebble velocity at r peak and the normalized di GLYPH<11> usivity at this location, GLYPH<15> = v peb( r peak) = D ( r peak) r peak. Initially, we start with an estimate for GLYPH<15> , which can, for example, be taken from the zero-model solution (Sect. 2.7). For a certain GLYPH<15> the semi-analytical model then updates the surface density at the peak ( GLYPH<6> peak) as well as r peak (or rather the di GLYPH<11> erence between r ice and r peak). With these updates we obtain a new value of GLYPH<15> and the procedure can be repeated until the fractional change in the parameters has become su GLYPH<14> ciently low.  \nWe emphasize that the semi-analytical model gives approximate solutions (correct at the GLYPH<24> 20% level), but that they are very useful since they come at almost zero computational costs. It is therefore ideal to be applied to the parameter searches that we consider in this section. In Appendix C the model is described in detail and compared to runs of the numerical steady-state model.", '5.2. Enhancement of the ice surface density': 'We first look at the relative e GLYPH<11> ect of enhancing the ice surface density near the snowline: for each disk model, we normalize the resulting peak ice surface density by the ice surface density at the peak radius in case there is only advection:  \nf GLYPH<6> ; peak = GLYPH<6> ice ; peak (1 GLYPH<0> GLYPH<16> ) F s = g ˙ M gas = 2 GLYPH<25> r peak v peb (35)  \nwhere v peb is calculated in the absence of condensation and evaporation. As before, we take GLYPH<16> = 0 : 5. In Fig. 7a we plot contours of constant f GLYPH<6> ; peak as function of GLYPH<11> and F s = g, where ˙ M gas is fixed at 10 GLYPH<0> 8 M GLYPH<12> yr GLYPH<0> 1 and GLYPH<28> 3 is fixed at 0 : 03. Clearly, the higher the solids-to-gas accretion rate F s = g, the higher f GLYPH<6> ; peak - at fixed GLYPH<11> the normalized ice peak increases from low F s = g-values to high F s = gvalues, due to collective e GLYPH<11> ects playing an increasingly important role. The plot also shows that disks with GLYPH<11> -values of about 10 GLYPH<0> 3 are best at enhancing the ice surface density with respect to the advection-only expected ice surface density. This is because higher values of GLYPH<11> lead to a broader and relatively lower peak,  \nFigure 7b shows contours of constant f GLYPH<6> ; peak, as function of GLYPH<11> and ˙ M gas. F s = g is fixed at 0.4 and GLYPH<28> 3 is fixed at 0.03. The cyan lines correspond to a Toomre parameter ( QT ) of unity at 10 au (1 au). QT is defined as:  \n<!-- image -->  \n<!-- image -->  \nFig. 6. Comparison of steady-state results for the simple, single-seed, fiducial model between the case where vapor has the same vertical distribution as the background gas ( HZ = H gas) and the case where vapor has the same scale height as the pebbles ( HZ = H peb). In all three panels, solid lines correspond to the HZ = H gas case and dashed lines correspond to the HZ = H peb case. (a) Steady-state water ice (blue) and water vapor (green) surface densities. (b) Midplane solids-to-gas ratios (green) and solids-to-gas surface densities ratios (blue). (c) Typical pebble mass m p.  \n<!-- image -->  \n<!-- image -->  \nF  \n<!-- image -->  \n/circledot  \nFig. 7. (a) Contours of constant relative ice surface density peak heights f GLYPH<6> ; peak (see text for details) as function of turbulence parameter GLYPH<11> and solids-to-gas accretion rate F s = g. The gas accretion rate ˙ M gas is fixed to 10 GLYPH<0> 8 M GLYPH<12> yr GLYPH<0> 1 and the initial size of the pebbles GLYPH<28> 3 is fixed to 0 : 03. (b) Contours of constant relative ice surface density peak heights (see text for details) as function of the turbulence parameter GLYPH<11> and the gas accretion rate ˙ M gas. The solids-to-gas accretion rate ratio F s = g is fixed to 0.4 and the initial size of the pebbles GLYPH<28> 3 is fixed to 0 : 03. The cyan lines denote a Toomre parameter ( QT ) of unity at 10 au (1 au); the space below these lines has QT < 1 and is gravitationally unstable at 10 au (1 au).  \nwhile lower GLYPH<11> -values, indicating larger gas densities, move pebbles into the Stokes drag regime. The fact that the relative height of the ice peak becomes smaller when pebbles enter the Stokes regime, is explained by the pebble size-dependency of the drift velocity. In the Epstein regime the drift velocity depends linearly on the pebble size, whereas in the Stokes regime the drift velocity is proportional to the square of the pebble size. This means that when pebbles enter the Stokes regime, their stopping time - and hence their drift velocity - increases rapidly. The e GLYPH<11> ect of this is that f GLYPH<6> ; peak decreases rapidly. Therefore, the transition from the Epstein regime to the Stokes drag regime is characterized by closely-spaced contours.  \nQT = cs GLYPH<10> GLYPH<25> G GLYPH<6> gas (36)  \nwhere G is the gravitational constant. A disk is gravitationally unstable if QT < 1. In Fig. 7b, models in the parameter space below the cyan line have QT < 1 at 10 au (1 au). In Fig. 7b  \ndensely-spaced contours again indicate the transition between Epstein and Stokes drag.  \nThe enhancement of the ice surface density near the snowline is relatively modest, but enrichment factors of a few can still be important for triggering streaming instability (Johansen et al. 2009). Besides, Fig. 7 corresponds to our single-seed model, which is conservative. Larger enhancements can be achieved in the framework of the many-seeds model (see Fig. 4). We can imitate the many-seeds model in the analytical model of Sect. C, crudely, when we let GLYPH<16> ! 0, which assumes that micron-sized silicate grains behave like vapor. In that case we see f GLYPH<6> ; peak increase by another factor of two (results not plotted). Therefore, total enhancements of f GLYPH<6> ; peak GLYPH<24> 5-10 are feasible, which significantly boost the likelihood of triggering streaming instability.', '5.3. Solids-to-gas ratios': "Let us now look at the absolute results, in terms of the midplane solids-to-gas ratio. In Fig. 8a, black lines denote contours of constant peak midplane solids-to-gas ratios, as function of GLYPH<11> and F s = g. We again fix ˙ M gas = 10 GLYPH<0> 8 M GLYPH<12> yr GLYPH<0> 1 and GLYPH<28> 3 = 0 : 03. The orange contours correspond to the total amount of solids in units  \nα  \nF  \nFigure 8b is the same as Fig. 7b, except that the black contours now denote constant values of the peak midplane solidsto-gas ratio. Just like in Fig. 8a, orange contours denote constant M solids ; needed, but now with ˙ M gas as the parameter on the x -axis. The region with closely-spaced contours in Fig. 8a and the oblique regions in the otherwise fairly horizontal contours in Fig. 8b are again due to the transition between the Epstein and Stokes drag regimes.  \n<!-- image -->  \n<!-- image -->  \n/circledot  \nFig. 8. (a) Contours of midplane solids-to-gas ratios (black) as function of the turbulence parameter GLYPH<11> and the solids-to-gas accretion rate F s = g. The gas accretion rate ˙ M gas is fixed to 10 GLYPH<0> 8 M GLYPH<12> yr GLYPH<0> 1 and the initial size of the pebbles GLYPH<28> 3 is fixed to 0 : 03. Orange contours denote the amount of solid material (in Earth masses) that is required to form the eventual peak in the ice surface density just outside the snowline. (b) Contours of midplane solids-to-gas ratios as function of GLYPH<11> and ˙ M gas. The solids-to-gas accretion rate F s = g is fixed to 0.4 and the initial size of the pebbles GLYPH<28> 3 is fixed to 0 : 03. Orange contours denote the amount of solid material (in Earth masses) that is required to form the eventual peak in the ice surface density just outside the snowline. The cyan lines denote a Toomre parameter ( QT ) of unity at 10 au (1 au); the space below these lines has QT < 1 and is gravitationally unstable at 10 au (1 au).  \nof Earth masses that is needed to form the eventual peak. This quantity ( M solids ; needed) is calculated as:  \nM solids ; needed GLYPH<24> F s = g ˙ M gas( r peak GLYPH<0> r snow) v gas (37)  \nwhere ( r peak GLYPH<0> r snow) =v gas is the typical timescale on which the peak forms, as discussed in Sect. 4.1. The value of M solids ; needed is similar to the amount of solids that gets 'locked up' in the peak 6 . Figure 8a shows that the higher the value of F s = g, the easier streaming instability can be triggered, as one may expect. It also shows that in order to reach midplane solids-to-gas ratios of GLYPH<24> 10%, a solids-to-gas accretion rate F s = g of the same order is required. This also holds for other values of the gas accretion rate ˙ M gas and initial pebble size GLYPH<28> 3.  \nUp to this point we have kept the parameter GLYPH<28> 3 constant at its fiducial value of 0.03. Recently, Yang et al. (2016) have shown that the streaming instability mechanism can also work for small particles ( GLYPH<28> GLYPH<24> 10 GLYPH<0> 3 ). In Fig. 9 we now vary GLYPH<28> 3 and plot the resulting peak midplane solids-to-gas ratio, whilst fixing ˙ M gas at 10 GLYPH<0> 8 M GLYPH<12> yr GLYPH<0> 1 , F s =g at 0.4, and GLYPH<11> at 3 GLYPH<2> 10 GLYPH<0> 3 . Di GLYPH<11> erent regimes in the plot are denoted by Roman numbers I-IV. In region I, the particles are well-coupled to the gas. The solids-to-gas ratio flattens out towards low GLYPH<28> 3-values, because radial drift and settling do not play a role for these particles. The boundary between regime I and regime II corresponds to the value of GLYPH<11> . The  \nhighest solids-to-gas ratios are reached for particles with stopping times of the same order as GLYPH<11> or smaller. For higher values of GLYPH<28> , the larger drift velocity (proportional to GLYPH<28> ) starts to dominate over the larger degree of settling (proportional to p GLYPH<28> ). Therefore, in regime II, the solids-to-gas ratio decreases with increasing GLYPH<28> . The transition from regime I and regime II corresponds to the GLYPH<28> 3-value for which GLYPH<15> is unity. The GLYPH<15> parameter (introduced in Sect. 5.1) increases with increasing GLYPH<28> 3, and the peak surface density decreases at a faster rate with increasing GLYPH<15> for GLYPH<15> > 1 than for GLYPH<15> < 1 (Eq. (C.13)). The physical explanation is that the solution changes from di GLYPH<11> usion-dominated ( GLYPH<15> < 1) to driftdominated ( GLYPH<15> > 1). Therefore, the solids-to-gas ratio decreases more steeply in regime III than in regime II. The final transition, from regime III to regime IV, is caused by the transition from the Epstein regime to the Stokes regime. For GLYPH<28> 3 values of GLYPH<24> 1, the enhancement e GLYPH<11> ect becomes minimal and the peak solids-to-gas ratio flattens out.  \nFrom Fig. 8a and Fig. 8b we learn that in the context of the GLYPH<11> -viscosity disk model, the largest solids-to-gas ratios near the snowline are reached for high GLYPH<11> -values. This means that our diffusion model does not require a very quiescent disk to achieve high solids-to-gas ratio, and also works in the early disk phase when GLYPH<11> and ˙ M gas are presumably high. Our results show that for high GLYPH<11> and high ˙ M gas, tens of Earth masses are locked up in pebbles in an annulus outside the snowline with a width of the order of 1 au.", '6.1. Conditions for streaming instability': 'Our results show that the ice surface density can be enhanced by a factor GLYPH<24> 3-5 outside the snowline, due to the e GLYPH<11> ect of water di GLYPH<11> usion and subsequent condensation onto icy pebbles (adopting GLYPH<16> < 0 : 5 or the many-seeds model increases the enhancement by another factor GLYPH<24> 2). This di GLYPH<11> ers from the results of Stevenson & Lunine (1988), who found that the same e GLYPH<11> ect could lead to an enhancement by as much as 75. However, they considered  \nα  \nFig. 9. Peak midplane solids-to-gas ratio as function of GLYPH<28> 3. ˙ M gas is fixed at 10 GLYPH<0> 8 M GLYPH<12> yr GLYPH<0> 1 , F s =g at 0.4, and the value of GLYPH<11> is 3 GLYPH<2> 10 GLYPH<0> 3 . We can distinguish four regimes I-IV, which are explained in the main text.  \n<!-- image -->  \na static system, whereas our model accounts for radial drift of pebbles and gas accretion onto the star. Even though we found a modest enhancement, we showed that the pebble density at the midplane can reach tens of percent of the gas density, under the condition that the solids-to-gas accretion rate F s = g is of order 0.1.  \nBut what are realistic values for F s = g? Ida et al. (2016) have calculated the pebble mass flux ˙ M peb in a viscously-evolving disk, taking into account dust growth and radial drift. They found that after 10 6 yr, F s = g GLYPH<24> 0 : 3 for a turbulence parameter GLYPH<11> = 10 GLYPH<0> 3 and a gas accretion rate ˙ M gas = 10 GLYPH<0> 8 M GLYPH<12> yr GLYPH<0> 1 . Since this result suggests F s = g-values of order 0.1 are plausible, water condensation can plausibly trigger streaming instability near the snowline. According to Ida et al. (2016), ˙ M peb is inversely proportional to GLYPH<11> . This implies that even though our results show that for high GLYPH<11> -values, water di GLYPH<11> usion and condensation leads to the largest solids-to-gas ratio at the snowline, a smaller pebble mass flux is expected, which reduces the e GLYPH<11> ect. The pebble mass flux is, however, dependent on time t as ˙ M peb / t GLYPH<0> 1 = 3 , meaning that the pebble mass flux earlier in the disk is higher than the fiducial 30% of the gas accretion rate quoted before. Because higher GLYPH<11> -values lead to larger solids-to-gas ratios but to lower pebble mass fluxes, we conclude that intermediate GLYPH<11> -values ( GLYPH<24> 10 GLYPH<0> 3 ) are most favorable to trigger streaming instabilities. The fact that the di GLYPH<11> usion-driven model of this work already produces enhancement of solids at intermediate GLYPH<11> -values and high ˙ M gas, implies that planetesimal formation does not have to await the arrival of quiescent disk conditions. Provided a su GLYPH<14> ciently large pebble flux, our results imply that planetesimals can form early.', '6.2. Interior or exterior?': 'In Saito & Sirono (2011), it is argued that the evaporation of icy pebbles leads to a large enhancement of small silicate grains interior to the snowline. In their model, micron-sized silicate grains are instantaneously released at the snowline, and then follow the radial motion of the gas. Since the gas moves at a much lower velocity than the drift velocity of the icy pebbles, the silicate grains pile-up and become gravitationally unstable. A similar mechanism for the formation of planetesimals in the inner disk has recently been presented by Ida & Guillot (2016).  \nWe confirm that the solids-to-gas ratio interior to the snowline is boosted due to the slow radial motion of the small silicates (see Fig. 4c), but find smaller enhancements than Saito & Sirono (2011) and Ida & Guillot (2016). These works assumed that the evaporated silicate grains remained in the same vertical layer as the icy pebbles, which had settled to the midplane due to their larger size. This is a natural assumption if one adopts instantaneous evaporation. However, in our model silicate grains are not instantaneously released from the pebbles at a single location, but are instead released across the ice evaporation front - the region between the snowline and the peak radius. With the exception of very low GLYPH<11> -values, the width of the ice peak is (much) larger than the gas scale height. Consequently, the timescale for vertical mixing is smaller than the timescale on which the silicate grains traverse the evaporation front, justifying our assumption that silicate grains (and water vapor) are distributed over H gas. Allowing for the vertical mixing of silicate grains is why we found smaller solids-to-gas ratios interior to the snowline than Saito & Sirono (2011) and Ida & Guillot (2016).  \nIn this work we did not model the porosity of pebbles, but assumed pebbles are compact spheres. For homogeneous porous spheres, the results are the same as for compact spheres with the same stopping time, if they are in the Epstein regime. This is because the evaporation and condensation rates (Eq. (13)-(14)) depend on s 2 p = m p (the surface to mass ratio of a pebble), which is e GLYPH<11> ectively the Epstein stopping time (Eq. (7)). For fractal aggregates, however, the evaporation rate could conceivably be increased due to a larger total surface area available for evaporation. In that case, the location of the snowline ( r snow) and of the peak ( r peak) are pushed to larger radial distances, because drifting pebbles start to evaporate earlier. The shape of the solids distribution profile and hence the width of the evaporation front, however, remain similar. We have confirmed this behaviour by running our simulation with an evaporation rate that has been arbitrarily increased by a factor 100 compared to Eq. (13). Therefore, a fractal interior structure of the pebbles would not much change the extent of the region in which pebbles evaporate. We note however that in the many-seeds model, a large increase in opacity due to small silicate grains might lead to a steeper temperature profile than assumed here, and therefore to a more narrow evaporation front.  \nEven though we found smaller enhancements interior to the snowline, adopting the many-seeds model leads to larger enhancements outside the snowline compared to the single-seed model. This is because outward radial di GLYPH<11> usion and subsequent sticking of silicate grains to icy pebbles adds to the ice peak and therefore enforces the enhancement of the solids-to-gas ratio outside the snowline (see Fig. 5).', '6.3. Condensation onto small grains': 'In our single-seed model, we adopted a single-size assumption for the solid particles, applicable to a drift-limited scenario, in which particles grow until they decouple from the gas and start to drift inward (see, e.g. , Birnstiel et al. (2010b); Krijt et al. (2016b)). Because the solids density is larger close to the star than in the outer disk, particles close to the star start to drift inward earlier than particles in the outer disk. This leads to an inside-out clearing of the solids component of the disk. In this case, one would not expect small particles to be present in significant amounts near the snowline at the time when icy pebbles approach the snowline, drifting in from the outer disk. Therefore, in the drift-limited case the single-size approximation for pebbles is valid.  \nThe single-size assumption does not hold for a fragmentation-limited size distribution, however. In that case, pebbles dominate the total mass of solids, but small particles dominate the total surface area (see, e.g. , Birnstiel et al. (2012)). Then, condensation of water vapor will predominantly happen onto small particles, instead of onto pebbles as considered in this work. Also, in our many-seeds model, we do not allow for condensation onto the small grains that get released from icy pebbles upon evaporation.  \nWithin the single-seed framework, we can mimic the situation where condensation happens onto small grains rather than onto pebbles by adopting a very small typical stopping time just outside the snowline ( GLYPH<28> 3). This results in approximately the same or even larger peaks, as illustrated by Fig. 9. In that sense, assuming condensation occurs solely onto pebbles, rather than onto small grains, is a conservative choice. However, a more correct approach would allow for condensation onto small grains in the many-seeds model, in which icy pebbles break up into small grains upon evaporation, or, alternatively, within a fragmentation-limited scenario by adopting a full size distribution. Such improvements will be considered in future work.', '6.4. Relevance of coagulation': "The many-seeds model ignores coagulation among the small silicate grains, as well as the possibility that water vapor condenses on to silicates. Because there are so many of them, the average silicate grain size would not increase significantly due to water condensation. In that case, the assumption that they are wellcoupled to the gas would still hold. However, growth of silicates might also occur due to coagulation. Growth increases the drift velocity, thereby reducing the e GLYPH<11> ect of the outward-di GLYPH<11> usion and recondensation of the silicate component. Consequently, coagulation among silicates would bring the many-seeds model closer to the single-seed model due to a decrease in the enhancement of 'locked' silicates outside the snowline. A similar trend is expected when the silicates that are released upon evaporation are larger. According to the experimental work of Aumatell & Wurm (2011), an icy aggregate of 1 cm in size breaks up into about 200 sub-aggregates when it evaporates. In any case, we expect that a more realistic 'break-up' or coagulation model will lead to results that lie in between the results for our many-seeds and single-seed model designs.  \nIn our work, we have also neglected coagulation among pebbles. On the one hand, one could think that icy pebbles on their way to the snowline would grow even more because of coagulation - at least as long their relative velocities stay low enough. On the other hand, according to Sirono (2011) and Okuzumi et al. (2016), the fragmentation threshold velocity for icy pebbles near the snowline is reduced due to sintering. This might motivate a smaller average icy pebble size (a lower GLYPH<28> 3 in our nomenclature) and therefore a smaller radial drift velocity, which would further increase the solid surface density exterior to the snowline (see Fig. 9).", '6.5. Observational implications': "We showed that intermediate-to-high GLYPH<11> -values are most favorable for triggering streaming instability near the snowline. The onset of streaming instability manifests itself through the formation of a peak in the solids-to-gas ratio that grows in height (increasing GLYPH<6> ice) but also in width - the latter because collective effects cause pebbles to move slower. Even though more solids are  \nrequired for high GLYPH<11> and high ˙ M gas than for low GLYPH<11> and low ˙ M gas, tens of Earth masses are su GLYPH<14> cient to form the ice peak. This also means that tens of Earth masses can be stored in an annulus outside the snowline, especially when the system is characterized by high pebble and gas accretion rates. Such annuli of solid enhancement can be observable with facilities as ALMA. In the context of the snowline, the ring structure seen in the TW Hya system with a width of about 1 au (Andrews et al. 2016) might very well correspond to a realisation of the water-di GLYPH<11> usion e GLYPH<11> ect in a highly turbulent disk with a high gas accretion rate.  \nThe fact that the planets in the inner Solar System are waterpoor (Marty 2012; McCubbin et al. 2012), even though the snowline should have migrated to 1 au before the disk was depleted, was pointed out as the 'snowline problem' by Oka et al. (2011). The early formation of a protoplanet near the snowline might provide a solution to this problem, as proposed recently by Morbidelli et al. (2016). Once a protoplanet near the snowline has formed it could halt the inward-drifting pebble flux, either by accreting the pebbles (Guillot et al. 2014), or by trapping them in pressure maxima created by the newborn planet (Zhu et al. 2014; Lambrechts et al. 2014). In the meantime, the water vapor gets accreted to the star and the snowline location migrates inward ( e.g. , Oka et al. (2011). Consequently, an early formation of a protoplanet near the snowline might explain the lack of water in the inner Solar System. The water di GLYPH<11> usion / condensation e GLYPH<11> ect discussed in this paper provides a way to form planetesimals near the snowline in an early stage of the disk ( i.e. , when the snowline is still located outside of the current water-poor region of the Solar System). After planetesimals of a certain size have formed, they can subsequently accrete pebbles and quickly grow to larger bodies ( e.g. , Ormel & Klahr (2010); Lambrechts & Johansen (2012, 2014); Bitsch et al. (2015); Visser & Ormel (2016)). Therefore, our water di GLYPH<11> usion / condensation model provides the first key step to realize the early formation of a protoplanet near the snowline.", '7. Conclusions': "Our main findings can be summarized as follows.  \n- 1. Water di GLYPH<11> usion and condensation near the snowline can result in an enhancement of a factor several in the ice surface density. Because of the dynamic setup of our model where water vapor is carried away with the accreting gas, this enhancement is much less than found by Stevenson & Lunine (1988). Nevertheless, the boost in the solids-to-gas ratio can still trigger streaming instability, provided a large enough pebble flux.\n- 2. The peak in the ice surface density is not located at the snowline, but exterior to it ( r peak > r snow). With larger incoming pebble flux, the peak becomes broader, because the backreaction of the solids on the gas reduces the pebble drift velocity. Depending on the turbulence strength, the width of the peak can be as much as GLYPH<24> 1 au.\n- 3. Such broad peaks can contain tens of Earth masses in pebbles, appearing as bright annuli at radio wavelengths.\n- 4. The release of many micron-sized silicate grains upon evaporation of icy pebbles produces a peak in 'locked' silicates exterior to the snowline, due to di GLYPH<11> usion and sweep-up of the silicates. This peak adds to the ice peak and therefore enforces the enhancement of the solids-to-gas ratio outside the snowline.  \n- 5. Interior to the snowline, the solids-to-gas ratio is boosted if many micron-sized silicate grains are released during evaporation, because they cause a 'tra GLYPH<14> c jam' e GLYPH<11> ect. However, because silicate grains mix with the background gas before they cross the snowline, the solids-to-gas ratio enhancement just interior to r snow is limited.\n- 6. In the context of a viscous disk model, the ratio between the accretion rate of solids and gas needs to be of order GLYPH<24> 0 : 1 in order to reach solids-to-gas midplane ratios of the order of tens percent. The mechanism operates best at intermediate ( GLYPH<24> 10 GLYPH<0> 3 ) GLYPH<11> -values. Therefore, planetesimals can form at an early time in the evolution of the disk.  \nAcknowledgements. D.S. and C.W.O are supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). We would like to thank Jonathan Guyer and Daniel Wheeler for their elaborate answers to FiPy -related questions. We would also like to thank Carsten Dominik, Shigeru Ida, Satoshi Okuzumi, Anders Johansen, Daniel Carrera and the other participants of the Lorentz workshop 'New directions in Planet Formation' (Leiden, July 2016) for stimulating discussions. We are indebted to Sebastiaan Krijt, Je GLYPH<11> Cuzzi, and the anonymous referee for constructive feedback on the manuscript.", 'References': 'ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3  \nAndrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40 Armitage, P. J. 2015, ArXiv e-prints [ arXiv:1509.06382 ] Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2 Aumatell, G. & Wurm, G. 2011, MNRAS, 418, L1 Bai, X.-N. & Stone, J. M. 2010, ApJ, 722, 1437 Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15 Benz, W. 2000, Space Sci. Rev., 92, 279 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010a, A&A, 513, A79 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010b, A&A, 513, A79 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112 Blum, J. & Wurm, G. 2000, Icarus, 143, 138 Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 Chapman, S. & Cowling, T. 1970, The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Di GLYPH<11> usion in Gases, Cambridge Mathematical Library (Cambridge University Press)  \nCiesla, F. J. & Cuzzi, J. N. 2006, Icarus, 181, 178  \nCuzzi, J. N. & Hogan, R. C. 2003, Icarus, 164, 127  \nCuzzi, J. N. & Zahnle, K. J. 2004, ApJ, 614, 490  \nDominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647  \nDra¸˙zkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105  \nDra¸˙zkowska, J. & Dullemond, C. P. 2014, A&A, 572, A78  \nEstrada, P. R. & Cuzzi, J. N. 2008, ApJ, 682, 515  \nEstrada, P. R., Cuzzi, J. N., & Morgan, D. A. 2016, ApJ, 818, 200  \nGuillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72  \nGuyer, J. E., Wheeler, D., & Warren, J. A. 2009, Computing in Science & Engineering, 11, 6  \nIda, S. & Guillot, T. 2016, A&A, 596, L3  \nIda, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72  \nJohansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022  \nJohansen, A. & Youdin, A. 2007, ApJ, 662, 627  \nJohansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75  \nKenyon, S. J. & Hartmann, L. 1987, ApJ, 323, 714  \nKretke, K. A. & Lin, D. N. C. 2007, ApJ, 664, L55  \nKrijt, S. & Ciesla, F. J. 2016, ApJ, 822, 111  \nKrijt, S., Ciesla, F. J., & Bergin, E. A. 2016a, ArXiv e-prints [ arXiv:1610.06463 ]  \nKrijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016b, A&A, 586, A20  \nLambrechts, M. & Johansen, A. 2012, A&A, 544, A32  \nLambrechts, M. & Johansen, A. 2014, A&A, 572, A107  \nLambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35  \nLichtenegger, H. I. M. & Komle, N. I. 1991, Icarus, 90, 319  \nLodders, K. 2003, ApJ, 591, 1220  \nLynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603  \nMarty, B. 2012, Earth and Planetary Science Letters, 313, 56  \nMcCubbin, F. M., Hauri, E. H., Elardo, S. M., et al. 2012, in Lunar and Planetary Science Conference, Vol. 43, Lunar and Planetary Science Conference, 1121  \nMorbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418  \nNakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375  \nNomura, H., Tsukagoshi, T., Kawabe, R., et al. 2016, ApJ, 819, L7  \n- Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141\n- Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82  \nOrmel, C. W. 2014, ApJ, 789, L18  \nOrmel, C. W. & Cuzzi, J. N. 2007, A&A, 466, 413  \nOrmel, C. W. & Klahr, H. H. 2010, A&A, 520, A43  \nOrmel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215  \nPollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62  \nRos, K. & Johansen, A. 2013, A&A, 552, A137  \nSafronov, V. S. 1969, Evolution of the Protoplanetary Cloud and Formation of Earth and the Planets, ed. V. S. Safronov (Moscow: Nauka. Transl. 1972 NASA Tech. F-677)  \nSaito, E. & Sirono, S.-i. 2011, ApJ, 728, 20  \nSato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15  \nShakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337  \nSirono, S.-i. 2011, ApJ, 733, L41  \nStevenson, D. J. & Lunine, J. I. 1988, Icarus, 75, 146  \nTanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414  \nTesti, L., Birnstiel, T., Ricci, L., et al. 2014, Protostars and Planets VI, 339  \nVisser, R. G. & Ormel, C. W. 2016, A&A, 586, A66  \nWeidenschilling, S. J. 1977, MNRAS, 180, 57  \nWhipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211  \nYang,  \nC.-C.,  \nJohansen,  \nA.,  \n& Carrera,  \nD.  \n2016,  \nArXiv  \ne-prints  \narXiv:1611.07014  \n[  \n]  \nYoudin, A. N. & Goodman, J. 2005, ApJ, 620, 459  \nYoudin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588  \nZhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-n. 2014, ApJ, 785, 122  \nZsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57  \nTable A.1. List of frequently-used symbols', 'Appendix A: List of symbols': 'A list of frequently-used symbols is given in Table A.1.', 'Appendix B: Expression Collective Effects': 'Expressions for the aerodynamic drift of single particles have been derived by Whipple (1972) and Weidenschilling (1977) for non-accreting disks. Later, Nakagawa et al. (1986) accounted for the backreaction forces that the particles collectively exert on the gas. The gas motion is then outwards while the total dust + gas mass flux becomes zero. These expressions can be generalized to account for a size distribution of particles (Tanaka et al. 2005; Estrada & Cuzzi 2008; Bai & Stone 2010).  \nIn (steady) accretion disks gas moves inwards at a rate that is set by the viscosity. This modifies the single particle drift expressions. In particular, small particles, which show negligible radial drift, nevertheless move inwards as they are carried by the gas. While for the single particle case the modification of the drift expressions for accreting disks is straightforward, expressions accounting for both collective e GLYPH<11> ects and viscosity have, to the best of our knowledge, not yet been presented. We derive those here.  \nThe equations of motions in a frame co-rotating with the local Keplerian period read:  \nD v Dt = GLYPH<0> v GLYPH<0> u t stop GLYPH<0> 2 GLYPH<10> K GLYPH<2> v + F Euler GLYPH<0> dust (B.1a)  \nD u Dt = GLYPH<26> p GLYPH<26>g v GLYPH<0> u t stop GLYPH<0> 2 GLYPH<10> K GLYPH<2> u + F Euler GLYPH<0> gas + F r P + F GLYPH<23> (B.1b)  \nwhere v is the particle velocity and u the gas velocity, D = Dt is the material derivative, the first term on the RHS is the drag force, the second term is the Coriolis force, and F Euler = GLYPH<0> ( d GLYPH<10> K = dt ) GLYPH<2> r is the Euler force that arises due to the radial motion. It therefore appears in the radial equations. In addition, the equation of motion for the gas includes a hydrostatic correction due to a radial  \npressure gradient ( F r P ) and viscous forces ( F GLYPH<23> ). Following convention, we write the pressure gradient in terms of a nondimensional GLYPH<17> :  \nF r P = 1 GLYPH<26>g dP dr GLYPH<17> 2 GLYPH<17>v K GLYPH<10> K e r (B.2)  \nand use the thin disk approximation (consider column densities GLYPH<6> gas instead of GLYPH<26> ) for the viscous force:  \nF GLYPH<23> = 1 GLYPH<6> gas r GLYPH<1> T (B.3)  \nwhere T is the viscous stress tensor. In cylindrical coordinates, the only relevant term is:  \nTr GLYPH<30> = GLYPH<23> GLYPH<6> gas r d GLYPH<10> K dr ; (B.4)  \nwith which the viscous force becomes  \nF GLYPH<23> = 1 GLYPH<6> gas r GLYPH<1> T = 1 GLYPH<6> gas 1 r 2 @ @ r r 2 GLYPH<23> GLYPH<6> gas r d GLYPH<10> K dr ! e GLYPH<30> (B.5)  \nwhere GLYPH<23> is the kinematic viscosity. In the GLYPH<11> -disk model we have that GLYPH<23> GLYPH<6> gas is constant and therefore  \nF GLYPH<23> = GLYPH<0> 3 GLYPH<23> GLYPH<10> K 4 r e GLYPH<30> : (B.6)  \nIn order to solve for the drift velocities, we put the accelerations on the LHS of Eq. (B.1) zero, D u = Dt = D v = Dt = 0. This is justified, because a change in the drift velocities occurs on a timescale of GLYPH<24> r =v r , which is always much longer than the stopping time. Inserting Eq. (B.2) and Eq. (B.6) into Eq. (B.1) gives a system of four equations and four unknowns (the velocities). Its solution is:  \n0 B B B B B B B B B B B @ v r vGLYPH<30> ur u GLYPH<30> 1 C C C C C C C C C C C A = 1 GLYPH<28> 2 s + (1 + GLYPH<24> ) 2 0 B B B B B B B B B B B B @ GLYPH<0> 2 GLYPH<17>GLYPH<28> s GLYPH<0> 3 2 ˜ GLYPH<23> (1 + GLYPH<24> ) 3 4 ˜ GLYPH<23>GLYPH<28> s GLYPH<0> GLYPH<17> (1 + GLYPH<24> ) 2 GLYPH<17>GLYPH<28> s GLYPH<24> GLYPH<0> 3 2 ˜ GLYPH<23> (1 + GLYPH<28> 2 s + GLYPH<24> ) GLYPH<0> 3 4 ˜ GLYPH<23>GLYPH<28> s GLYPH<24> GLYPH<0> GLYPH<17> (1 + GLYPH<28> 2 s + GLYPH<24> ) 1 C C C C C C C C C C C C A v K (B.7)  \nwhere GLYPH<28> s = t stop GLYPH<10> K , GLYPH<24> = GLYPH<26> d =GLYPH<26>g , and ˜ GLYPH<23> = GLYPH<23>= r 2 GLYPH<10> K . Note that in the GLYPH<11> -prescription ˜ GLYPH<23> GLYPH<24> GLYPH<11>GLYPH<17> , that is, for GLYPH<11> GLYPH<28> 1 the modification due to viscosity is minor.', 'Appendix C: Semi-analytical model': 'In this section we aim for an analytic prescription of the key characteristics of the steady-state ice and vapor distribution for any given disc model - not necessarily with the same gas and temperature profiles as assumed in this study. We will derive expressions for the peak of the ice surface density and its location, which enable us to determine whether the conditions for streaming instability beyond the snowline can be reached.', 'Appendix C.1: Location of the snowline': 'The location of the snowline, r snow, is the outer-most radius where the entire incoming ice flux can thermodynamically exist in the gaseous state. Its location is found by equating the equilibrium (saturation) pressure to the vapor pressure of H2O corresponding to the incoming icy pebble flux, that is, P eq = PZ ; a,  \nwhere PZ ; a is the steady-state vapor pressure; i.e. , the vapor pressure obtained when all the H2O is in the gas phase. Therefore, the snowline r snow is given by the radius r where  \nP eq ; 0 exp " GLYPH<0> Ta T ( r ) # = kBT GLYPH<22> GLYPH<6> Z ; a p 2 GLYPH<25> H gas = kBT ( r ) GLYPH<22> ˙ M ice = 3 GLYPH<25>GLYPH<23> ( r ) p 2 GLYPH<25> H gas( r ) (C.1)  \nsee Eq. (11) and Eq. (34). Given ˙ M ice, GLYPH<11> and a disk model, we can therefore readily solve (numerically) for r snow.  \nA caveat of Eq. (C.1) is that the role of transport processes are neglected. These could potentially introduce disequilibrium corrections. However, we find that the times to achieve the equilibrium (saturation) pressure is typically much shorter than the transport timescales (by di GLYPH<11> usion or drift).', 'Appendix C.2: The equilibrium density': 'Generalizing the above arguments to locations beyond r snow, we approximate the vapor pressure with the saturation pressure - an approximation that is again justified as long as the timescales for evaporation (and condensation) are su GLYPH<14> ciently short (in other words: the super-saturation level is small). Similarly, we can define a vapor surface density ( GLYPH<6> eq) corresponding to P eq. Rearranging Eq. (C.1) in terms of GLYPH<6> we obtain this equilibrium surface density:  \nGLYPH<6> eq GLYPH<17> p 2 GLYPH<25> H gas GLYPH<22> kBT P eq = GLYPH<6> snow r r snow ! p exp " Ta T snow 1 GLYPH<0> r r snow ! q !# (C.2)  \nwhere GLYPH<6> snow and T snow are the vapor surface density and the temperature at the location of the snowline. In Eq. (C.2) we have assumed power-laws for the temperature, gas surface density, and the gas scaleheight. In our case we have p = 7 = 4 and q = 1 = 2 (see Sect. 2.1).  \nTaylor-expanding the (weakly-varying) power-laws involved in the expression of GLYPH<6> eq around r snow gives: 7  \nGLYPH<6> eq( x ) GLYPH<25> GLYPH<6> snow (1 + px ) exp " GLYPH<0> qTax T snow # = GLYPH<6> snow (1 + px ) e GLYPH<0> a eq x (C.3)  \nwhere x = ( r GLYPH<0> r snow) = r snow is the fractional distance beyond the snowline and  \na eq = qTa T snow GLYPH<25> 15 GLYPH<18> q 0 : 5 GLYPH<19> GLYPH<18> T snow 200 K GLYPH<19> GLYPH<0> 1 (C.4)  \nis a dimensionless number that weakly depends on the location of the snowline.  \nIn Fig. C.1 the respective vapor density profiles have been plotted. Parameters correspond to those of the default model (Table 1). The simulated result of Sect. 4 is plotted by the thick gray curve. It is characterized by a sharp dent at r snow. Exterior to the snowline GLYPH<6> Z follows the equilibrium profile (solid black curve), while interior to r snow, it is limited by the imposed mass flux (solid dashed curve). The approximate expressions derived in Eq. (C.3), valid for r > r snow, are shown by the blue curves. We obtain r snow = 2 : 09 au, GLYPH<6> snow = 92 g cm GLYPH<0> 2 and a eq = 16 : 9. The solid blue curve approximates GLYPH<6> eq very closely. The dashed  \nFig. C.1. Vapor surface density profiles. Shown in black are the steady state vapor profile when all the ice is in the gas phase (thin solid curve), and the surface density corresponding to equilibrium vapor pressure P eq (Eq. (C.2); dashed solid curve). The snowline r snow is located at the intersection of these two curves. The numerically-obtained profile (thick gray curve) closely follows the minimum of steady-state and equilibrium profiles. The blue curves present mathematically-convenient approximations to GLYPH<6> eq (see text), valid for r > r snow.  \n<!-- image -->  \nblue curves gives a further approximation to GLYPH<6> eq, obtained by putting p = 0. This slightly underestimates GLYPH<6> Z but nevertheless gives a reasonable approximation that we will use below.  \nM Z GLYPH<25> M eq = GLYPH<0> v gas GLYPH<6> eq GLYPH<0> D gas GLYPH<6> gas d dr GLYPH<6> eq GLYPH<6> gas ! GLYPH<25> GLYPH<0> v gas GLYPH<6> eq GLYPH<0> D gas d GLYPH<6> eq dr (C.5)  \nAssuming GLYPH<6> Z = GLYPH<6> eq for the vapor density beyond r snow, we write for the associated mass flux M Z :  \nwhere we neglected the curvature of the disk, i.e., we assumed that the scales over which GLYPH<6> eq and GLYPH<6> gas change are much smaller than r . This is justified as long as x GLYPH<28> 1. In a similar vein, we assume that quantites as D gas and v gas do not vary. With Eq. (C.3) for GLYPH<6> eq we obtain:  \nM eq D gas = r snow = GLYPH<0> v gas r snow D gas GLYPH<6> eq GLYPH<0> d GLYPH<6> eq dx (C.6) = GLYPH<6> snow GLYPH<16> GLYPH<0> b gas(1 + px ) + a eq(1 + px ) GLYPH<0> p GLYPH<17> e GLYPH<0> a eq x GLYPH<25> GLYPH<6> snow( a eq GLYPH<0> b gas) e GLYPH<0> a eq x (C.7)  \nwhere b gas = v gas r snow = D gas is the dimensionless gas velocity and in the last step p = 0. In a viscously-evolving disk v gas = 3 GLYPH<23>= 2 r and hence b gas = 1 : 5.', 'Appendix C.3: The ice profile': 'With the equilibrium profile as the approximation to the vapor density, we solve for the surface density of the ice by invoking conservation of mass. In steady-state the total H2O mass flux M tot ; ice = M ice + M Z is constant. Interior to the snowline, M Z = M tot ; ice while for r GLYPH<29> r snow M Z GLYPH<28> M ice GLYPH<25> M tot ; ice. But near the snowline the mass fluxes change rapidly. Here, the steep gradient in the vapor density causes an outward (positive) mass flux, which increases the inwardly-directed icy pebble flux:  \nGLYPH<0>M ice GLYPH<17> GLYPH<6> ice v ice + D p d GLYPH<6> ice dr = M Z GLYPH<0> M tot ; ice : (C.8)  \nFig. C.2. Steady state surface density of ice for the default model parameters. The simulated profile from the numerical model (thick grey line) is well reproduced by the full analytical solution using p = 7 = 4 (solid curve; formula not shown in the main text) and the further p = 0 approximation (dashed curve; Eq. (C.10)). The black square denotes the values corresponding to the mass peak ( r peak ; GLYPH<6> peak) given by Eqs. (C.12) and (C.13).  \n<!-- image -->  \n(Note that just outside r snow, M ice and M tot ; ice are negative and M Z is positive.) With M eq for M Z Eq. (C.8), in terms of x , reads  \nd GLYPH<6> ice dx + b peb GLYPH<6> ice = r snow D p GLYPH<16> M eq GLYPH<0> M tot GLYPH<17> (C.9)  \nwhere b peb = r snow v peb = D p is the ratio between the speed at which the pebbles drift in and the radial velocity of the gas. Usually, pebbles outpace the gas and b peb GLYPH<29> 1.  \nTo solve Eq. (C.9) we will assume that b peb and other parameters are constant in x . Assuming constant b peb is clearly an approximation, since the pebbles may (i) acquire thick icy mantles on their approach to the snowline but (ii) lose all their ice once their are very close to r snow. Hence, the aerodynamical properties of the pebbles and therefore v peb are expected to vary considerably. Disregarding these concerns, momentarily, the solution to Eq. (C.9) reads:  \nGLYPH<6> ice = e GLYPH<0> b peb x Z x 0 e b peb x 0 r snow M eq( x 0 ) D p dx 0 = GLYPH<6> snow(1 GLYPH<0> GLYPH<15> gas) " c M (1 GLYPH<0> e GLYPH<0> a eq " x ) " + e GLYPH<0> a eq x GLYPH<0> e GLYPH<0> a eq " x " GLYPH<0> 1 # (C.10)  \nwhere we used Eq. (C.6) for M eq, put p = 0 8 and defined " = b peb = a eq, " gas = b gas = a eq and c M = GLYPH<0>M tot ; ice r snow = GLYPH<6> snow D p( a eq GLYPH<0> b gas) to keep the notation concise. Typically, " gas GLYPH<25> 1 : 5 = a eq and c M are small numbers. In a steady-state viscous disk GLYPH<6> snow GLYPH<25> ˙ M ice = 2 GLYPH<25> r v gas and M tot ; ice = GLYPH<0> ˙ M ice = 2 GLYPH<25> r snow. Hence in a viscous disk c M is related to " gas:  \nc M = v gas r snow D p( a eq GLYPH<0> b gas) = " gas 1 GLYPH<0> " gas GLYPH<25> 0 : 1 (C.11)  \nIn Fig. C.2 we plot GLYPH<6> ice for the parameters of the default model (Table 1). The simulated data are given by the thick grey curve while the analytical profiles are in blue - the more precise approximation (with p = 7 = 4; solid) and the p = 0 case in  \nTable C.1. Comparison between analytical and numerical model  \nF  \nNotes. Column entries denote: (S / F) simple (no GLYPH<22> -variations or collective e GLYPH<11> ects) or full model; (parameter) input parameter that is varied and its value; ( FiPy ) values for the position of the ice peak (in au) and the surface density (in cgs units) for the ice peak; (Analytical) same for the analytical model.  \ndashed. The analytical profiles fit the simulation data very well. The key parameter for the analytical profiles is the velocity of the pebbles v peb. Here, we have chosen v peb to be the pebble velocity at the ice peak r peak, v peb GLYPH<25> 2 : 4 m s GLYPH<0> 1 , which provides a decent fit. The fit is somewhat sensitive to the choice of v peb; for example, adopting v peb( r peak) with the initial value for m ice ( i.e. , without deposition of vapor) would overestimate the numerical curve by 20%.', 'Appendix C.4: Peak values and final tuning': 'Equation (C.10) has a maximum at ( x ; GLYPH<6> ice) = ( x peak ; GLYPH<6> peak) where  \nx peak = log [ " + c M ( " GLYPH<0> 1)] a eq( " GLYPH<0> 1) (C.12)  \nor r peak = (1 + x peak) r snow in dimensional units, and  \nGLYPH<6> peak = (1 GLYPH<0> " gas) c M + [ c M ( " GLYPH<0> 1) + " ] 1 = (1 GLYPH<0> " ) " GLYPH<6> snow : (C.13)  \nSince c M and " gas are fixed, the latter expression only depends on " . Lower " - meaning: a smaller pebble velocity at r peak - increases both the width ( x peak) and the magnitude of the ice peak, resulting in a more pronounced, \'fatter\' ice bump.  \nAll that remains is to find an expression for the pebble velocity at the ice peak v peb, which in turn depends on the amount of ice those pebbles have accreted. To obtain m peb = m core + m ice we invoke conservation of the pebble number flux  \nN p = ˙ M sil = m core 2 GLYPH<25> r v peb = GLYPH<6> ice m ice : (C.14)  \nHence, we can obtain m ice from GLYPH<6> peak, calculate the the aerodynamical properties of the pebbles at the peak ( i.e. , their stopping time) and obtain the pebble velocity v peb.  \nWith Eq. (C.14) we have a closed system of expressions to obtain the relevant quantities at the ice peak. These can best be computed by an iterative scheme, e.g. ,  \n- 1. From (an initial guess for) " obtain GLYPH<6> peak and the position of the ice peak r peak from the ice flux conservation model, Eq. (C.13).\n- 3. From m ice and the gas properties at r = r peak, compute the stopping time of the pebbles at r peak. Update the pebble velocity v peb and its normalized variant " .\n- 2. From GLYPH<6> peak, obtain the ice mass m ice at the ice peak from the pebble conservation law, Eq. (C.14).  \nThis scheme can easily be extended with collective and mean molecular weight e GLYPH<11> ects ( i.e. , the complete model). Finally, in order to obtain a better match to our numerical results we adopt two ad-hoc empirical \'fixes\'. First we reduce Eq. (C.13) by 10%, which accounts for the fact that the p = 0 approximation tends to overestimate the peak. Secondly, we slightly increase " when it drops below unity. This fix approximately accounts for curvature e GLYPH<11> ects (not included in the model) that become apparent when the ice peak becomes broad (at low " ). Hence we take " = max( " GLYPH<3> ; " GLYPH<3> 0 : 8 ) where GLYPH<15> GLYPH<3> = ( v peb ; peak r peak = a 2 D peak) is the uncorrected value.  \nIn Table C.1 a comparison of the model with the numerically-obtained peak parameters is given for a number of model parameters. Generally, the agreement is very good; the error in GLYPH<6> peak is at most 20%. Given its crudeness, the analytical model, however, does not always produce a perfect match - especially not when it comes to the profile. For example the low turbulence runs ( GLYPH<11> = 3 GLYPH<2> 10 GLYPH<0> 4 ), which imply high gas densities when ˙ M gas is kept the same, still show a good match to GLYPH<6> peak, but the profile corresponding to Eq. (C.10) is very di GLYPH<11> erent from the numerical model. The reason for this is that pebbles enter the Stokes drag regime when approaching the ice peak, causing a sharp increase in the stopping time (and " ). The stopping time at the ice peak is in that case not representative for other radii.'}