© The Authors 2025

1 Introduction

Planets are born in protoplanetary disks, going through a series of growth from dust grains to planetesimals prior to the formation of planet embryos and the final assembly of planets (Drążkowska et al. 2023; Birnstiel 2024). Tiny micron-sized (µm) grains, inherited from the molecular cloud, can readily grow via coagulation (Dominik & Tielens 1997; Ormel et al. 2007). However, when reaching centimeter-sized (cm) pebbles, growth is believed to stall due to two key factors. First, as relative velocities increase with dust size, collisions result in fragmentation or bouncing, rather than coagulation (Blum & Wurm 2000, 2008). Second, pebbles experience rapid inward radial drift, which quickly drains the solid mass budget (Weidenschilling 1977). To overcome these growth barriers, a local enhancement in solids is crucial. An elevated solid-togas ratio can trigger streaming instabilities (SI), which, in turn, induce the gravitational collapse of pebbles, leading to direct formation of km-sized planetesimals (Youdin & Goodman 2005; Johansen et al. 2007, 2009).

Disk snowlines, particularly the water snowline, are thought to be promising sites for planetesimal formation due to their capacity to enhance solid material (Drążkowska et al. 2023). These enhancements occur as a result of two effects. First, as pebbles drift across the water snowline, their ice contents sublimates, leaving behind pure silicates. This process creates a “traffic jam” effect if the decreasing stickiness of dry silicates results in a low fragmentation threshold (Dominik & Tielens 1997; Wada et al. 2013; Gundlach & Blum 2015) or if icerich pebbles disaggregate upon sublimation (Saito & Sirono 2011; Aumatell & Wurm 2011). Both scenarios would reduce the pebble size and slow down their inward radial drift, leading to a pile-up of silicates-dominated pebbles interior to the snowline (Saito & Sirono 2011; Hyodo et al. 2019, 2021). Second, exterior to the snowline, the density in ice can be boosted when vapor diffuses back across it and recondenses onto ice-rich pebbles, creating an additional pile-up (Cuzzi & Zahnle 2004; Ros & Johansen 2013; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). In both the “traffic jam” and “vapor retro-diffusion” mechanisms, numerical studies have shown that the midplane dust-to-gas ratio can readily reach unity (e.g., Schoonenberg & Ormel 2017; Hyodo et al. 2019), a level necessary to trigger SI (Johansen & Youdin 2007; Carrera et al. 2015; Bai & Stone 2010a; Li & Youdin 2021; Lim et al. 2024).

However, recently, the viability of the traffic jam effect has been called into question. Laboratory experiments find that ice’s stickiness drops sharply below 200 K, suggesting that silicates are not more fragile than ice (Gundlach et al. 2018; Musiolik & Wurm 2019). Furthermore, multi-band analysis of observations from the Atacama Large Millimeter/submillimeter Array (ALMA) towards V883 Ori, an outburst system that has extended its snowline to ≈40–80 au, reveals continuous dust sizes across the snowline (Houge et al. 2024). This leaves the vapor retro-diffusion mechanism as the key snowline process to elevate the solid-to-gas ratio.

Previous studies investigating factors conducive to solid pile-up typically employ 1D, vertically-averaged and isothermal assumptions. Under these assumptions, Schoonenberg & Ormel (2017) find that a high pebble flux, low Stokes number and large diffusivity are beneficial to form ice pile-up. Additionally, accounting for the back reaction of solids on the gas can double the enhancement of ice. The “traffic jam” of silicates can even act in a runaway fashion when considering the back reaction on diffusivity of dust and gas, which decreases with high solid concentration (Hyodo et al. 2019).

However, the vertical structure in the disk in reality plays an important role. First, in addition to radial transport, pebbles settle vertically and vapor can diffuse along the z-direction. Thus, recondensation-driven pebble growth could strengthen vertical settling, potentially enhancing pile-ups. Second, a vertical temperature gradient in the disk will determine the morphology of the snowline. The water snowline in the solar nebula might have resided in the actively-heated regions (e.g., Oka et al. 2011, but also see Mori et al. 2021), where the midplane is hotter than the upper layer due to viscous dissipation. In contrast, passively- heated disks usually hold cooler midplanes (e.g., Chiang & Goldreich 1997). These different vertical temperature structures lead to distinct snowline morphology as illustrated in Oka et al. (2011): active heating sublimates ice at the midplane, creating a narrow ice-condensing region above the midplane. Ros & Johansen (2013) demonstrated that a narrower ice-condensing region speeds up pebble growth by concentrating the material supply for recondensation. However, their Monte Carlo simulations, which were conducted in a local box with periodic boundaries in the radial direction, cannot directly address the consequences of snowline morphology to the solid pile-up.

Moreover, ice sublimation and vapor condensation involve latent heat energy exchange with the surrounding environment (e.g., Owen 2020). Although this is typically neglected in disk snowline modeling, latent heat absorption has been shown to significantly shift the sublimation front (i.e., the snowline), creating an “isothermal” region in accreting planets’ envelopes (Wang et al. 2023). Drążkowska & Alibert (2017) shows that different disk temperature profiles lead to different strength of pile-up and a stronger pile-up occurs when snowline is closer to the star due to the higher surface density there. Therefore, a 2D study that self-consistently includes pebble dynamics, vapor transport and phase change processes, while solving for the disk temperature and density structure, is required to assess the impact of the snowline on the characteristics of the solid pile-up.

A proper physical description of the snowline is also important for the post-planetesimal phases of planet formation. After planetesimals are formed in the pile-up, a substantial enhancement in solids would enable fast growth of planets (Ormel et al. 2017; Liu et al. 2019; Schoonenberg et al. 2019). Ormel et al. (2017) proposed that the snowline could act as a planetary “factory,” producing planets sequentially to form compact planetary systems such as Trappist-1 (Gillon et al. 2016, 2017). More generally, Jiang et al. (2023) demonstrates that dust rings, where solids are trapped, can spawn planets even at tens of au, where the ubiquitous rings are usually observed by ALMA (e.g., Huang et al. 2018). As the snowline delineates the solid and gas composition of the disk (e.g., Öberg et al. 2023), the density distributions of ice and vapor across it (e.g., the width, strength, and composition of the pile-up) will determine the composition of planets that form at the snowline region (Schoonenberg et al. 2019; Johansen et al. 2021).

In summary, given all these important roles of snowline in planet formation, a realistic, 2D description of the snowline structure will enhance our understanding of how the snowline influences the formation of the first generation of planetsi- mals, the growth of planetary embryos into planets, and the composition of the planets that ultimately emerge. Besides, self- consistently solving the temperature and density structures could help identify potential observation imprints of water snowlines, which still remain elusive (e.g., Zhang 2024).

In this work, we perform 2D multifluid hydrodynamic simulations in the disk’s radial (r)-vertical (z) plane with the inclusion of sublimation and condensation processes. Pebbles are modeled as compounds containing various chemical components (e.g., water ice, silicates). Each component is treated as a single pressureless dust fluid, whose dynamics is tracked using the multifluid dust module in Athena++ (Huang & Bai 2022). In addition, vapor species are treated as passive scalars. The mass transfer and latent heat energy exchange during sublimation or condensation are self-consistently solved with the recently- developed phase change module (Wang et al. 2023). This novel setup enables simultaneous tracking of gas flow, pebble dynamics and the phase change between ice and vapor in hydrodynamic simulations. While we have omitted, for simplicity, dust collisional processes (coagulation or fragmentation), the pebble size can be altered by vapor recondensation and ice sublimation. Furthermore, to determine the temperature structure near the snowline, we adopt a two-stream radiation transfer method (e.g., Mori et al. 2019) that takes into account stellar irradiation, viscous heating, heat diffusion, and latent heat exchange.

The paper is structured as follows. In Sect. 2, we present the 2D snowline model, covering the multifluid hydrodynamic equations for gas, vapor, and pebbles, phase change processes, and the radiation transfer method. Section 3 begins with describing the typical 2D steady-state snowline structure for the active, passive, and vertically-isothermal disk setups. Then, we present the impact of the snowline morphology on the solid pile-up, highlighting a “water-cycle” that promotes solid pile-up in active disks. In Sect. 4, we examine the effect of latent heat, which varies in different types of snowline thermal structure. Section 5 presents a detailed flow pattern analysis and Lagrangian water parcel tracking to understand the build-up of pile-up and the “water-cycle” phenomenon in 2D. In Sect. 6, we discuss the implications on planetesimal formation, pebble accretion, and the observational imprints left by latent heat effect. We further discuss the caveats and future improvements of this work. Section 7 gives our main conclusions.

2 Methods

This section describes the 2D hydrodynamic snowline model and presents the choices of simulation parameters. Section 2.1 outlines the gas disk model. Then, Sect. 2.2 presents the governing equations of the gas, pebble and vapor fluids, and the aerodynamic properties of pebbles are detailed in Sect. 2.3. The treatments of ice sublimation and vapor condensation are described in Sect. 2.4 and the radiation transfer method in Sect. 2.5. Finally, Sect. 2.6 presents the design of simulation cases explored in this work and summarize their corresponding parameter space.

2.1 Gas disk model

We set up a 2D disk model in the radial (r) and vertical (z) direction in a Cylindrical coordinate. Throughout the paper, the α-disk model (Shakura & Sunyaev 1973) is applied, where the viscosity is prescribed as $$\nu =\alpha c_s{H}_{\text{g}},$$(1)

where $c_s=\sqrt{k_{\text{B}}T/\mu m_{\text{p}}}$ is the local sound speed, µ is the mean molecular weight, H_g = c_s/Ω is the gas scale height and α is the nondimensional viscosity parameter. Assuming that the gas gets accreted to the central star at a steady rate Ṁ_acc , the gas surface density can be obtained as (e.g., Armitage 2020) $${\mathrm{\Sigma }}_{\text{g}}=\frac{{\dot{M}}_{\text{acc}}}{3\pi \nu }.$$(2)

In reality, the viscosity varies with temperature and therefore altitude in the disk. However, to simplify the disk profile, we fix the viscosity to be a power-law function of radius only, ν = ν₀(r/r₀)^b, where ν₀ = αc_s,0H_g,0 is the viscosity at the reference position r₀. We choose r₀ = 3.0 au and the power law index b = −1 following Schoonenberg & Ormel (2017). In this way, given the values of Ṁ_acc and α, the surface density profile can be fixed. The turbulent viscosity also causes diffusion to the gaseous components. We set the turbulent diffusivity, D_g , of vapor equal to the viscosity, ν.

2.2 Governing equations

The hydrodynamic equations of gas, pebble and vapor tracer fluids are presented in conservative form. The subscripts “g” and “p” are used to denote gas and pebble, respectively. We refer to Wang et al. (2023) for detailed design of this hydyodynamic system. The governing equations are expressed as: $$\frac{\partial {\rho }_{\text{g}}}{\partial t}+\nabla ·\left({\rho }_{\text{g}}{\nu }_{\text{g}}\right)=0,$$(3) $$\frac{\partial \left({\rho }_{\text{g}}{\nu }_{\text{g}}\right)}{\partial t}+\nabla \cdot \left({\rho }_{\text{g}}{\nu }_{\text{g}}{\nu }_{\text{g}}+P_{\text{g}}I+{\Pi }_{\nu }\right)={\rho }_{\text{g}}{f}_{\text{g},\text{src}},$$(4) $$\frac{\partial E_{\text{g}}}{\partial t}+\nabla \cdot \left[\left(E_{\text{g}}+P_{\text{g}}\right)v_{\text{g}}+{\Pi }_v\cdot v_{\text{g}}\right]={\rho }_{\text{g}}{f}_{\text{g},\text{src}}\cdot v_{\text{g}},$$(5) $$\frac{\partial {\rho }_{\text{p},i}}{\partial t}+\nabla \cdot \left({\rho }_{\text{p},i}{v}_{\text{p}}+{ℱ}_{\text{p},\text{dif},i}\right)=0,$$(6) $$\begin{array}{c}\frac{\partial {\rho }_{\text{p},i}\left(v_{\text{p}}+v_{\text{p},\text{dif},i}\right)}{\partial t}+\nabla \cdot \left({\rho }_{\text{p},i}{v}_{\text{p}}{v}_{\text{p}}+{\Pi }_{\text{dif},i}\right)\\ ={\rho }_{\text{p},i}{f}_{\text{p},\text{src},i}+{\rho }_{\text{p},i}\frac{v_{\text{g}}-v_{\text{p}}}{t_{\text{s}}},\hfill \end{array}$$(7) $$\frac{\partial {\rho }_{\text{vap}}}{\partial t}+\nabla \cdot \left({\rho }_{\text{vap}}{v}_{\text{g}}+{ℱ}_{\text{dif},\text{\hspace{0.17em}vap\hspace{0.17em}}}\right)=0.$$(8)

Here, v_g and v_p represent the velocities of gas and pebbles, respectively. Tracer fluid (passive scalar) has the same velocity as gas. The tensors Π_ν and Π_dif,i denote the viscous stress of gas and diffusion of momentum flux, respectively. For each component i (ice or silicate) of the pebble, the particle concentration flux is given by $${ℱ}_{\text{p},\text{dif},i}=-{\rho }_{\text{g}}{D}_{\text{p}}\nabla \left(\frac{{\rho }_{\text{p},i}}{{\rho }_{\text{g}}}\right)\equiv {\rho }_{\text{p},i}{v}_{\text{p},\text{dif},i},$$(9)

which defines the effective diffusive velocity v_p,dif,i. The last term in Equation (7) represents the aerodynamic drag experienced by the particles, which is controlled by the stopping time, t_s . The backreaction of solids acting on gas is neglected in this work (see Sect. 6.1). The phase change module tracks the vapor fraction using passive scalars, as shown in Equation (8). Finally, f_src represents the source term due to external forces, which only includes stellar gravity here. The hydrodynamic equations are solved in a spherical-polar coordinate with Athena++ (Stone et al. 2020), where the pebble dynamics is computed by the multifluid dust module (Huang & Bai 2022).

2.3 Pebble dynamics

Due to the aerodynamic drag of the gaseous disk, pebbles losing angular momentum go through inward radial drift. The aerodynamic properties of pebbles can be characterized by the stopping time (t_s). To obtain t_s , we account for two drag regimes: Epstein drag and Stokes drag. Then the stopping time is given by $$t_{\text{s}}=\{\begin{array}{l}\frac{{\rho }_{•,\text{p}}{s}_{\text{p}}}{v_{\text{th}}{\rho }_{\text{g}}},\left(\text{\hspace{0.17em}Epstein\hspace{0.17em}}:s_{\text{p}}<\frac{9}{4}{l}_{\text{mfp}}\right),\\ \frac{4{\rho }_{•,\text{p}}{s}_{\text{p}}^2}{9v_{\text{th}}{\rho }_{\text{g}}{l}_{\text{mfp}}},\left(\text{\hspace{0.17em}Stokes\hspace{0.17em}:\hspace{0.17em}\hspace{0.17em}}{s}_{\text{p}}>\frac{9}{4}{l}_{\text{mfp}}\right),\end{array}$$(10)

where the mean free path of the gas molecules l_mfp, $$l_{\text{mfp}}=\frac{\mu m_{\text{p}}}{\sqrt{2}{\rho }_{\text{g}}{\sigma }_{\text{mol}}},$$(11)

with σ_mol is the molecular collision cross section, which we approximate by the collision cross section of molecular hydrogen σ_mol = 2 × 10⁻¹⁵ cm² (Chapman & Cowling 1991), s_p is the pebble size, ρ_•,p is the pebble internal density and $v_{\text{th}}=\sqrt{8/\pi }{c}_{\text{s}}$ is the thermal velocity of the gas molecules.

We assume that pebbles drift towards the water snowline with a steady mass flux and consist of two material components: refractory silicate and volatile water ice. The water ice mass fraction of the particles, ζ = ρ_p,ice /(ρ_p,ice + ρ_p,sil) is assumed ζ = 0.5 for pebbles that enter the simulation domain. Therefore, pebbles are injected at a rate following Ṁ_peb = 2f_i/gṀ_acc, where f_i/g is the ice-to-gas flux ratio.

However, ζ is allowed to vary in the simulation due to sublimation and condensation processes. The internal density ρ_•,p of a pebble compound changes accordingly, $${\rho }_{•,\text{p}}={\left(\frac{1-\zeta }{{\rho }_{•}\text{\hspace{0.17em}sil\hspace{0.17em}}}+\frac{\zeta }{{\rho }_{•,\text{\hspace{0.17em}ice\hspace{0.17em}}}}\right)}^{-1},$$(12)

where ρ_•,p is the internal density of pebble. We adopt ρ_•,sil = 3.0 g cm⁻³ for pure silicate and ρ_•,ice = 1.0 g cm⁻³ for pure water ice.

To calculate the stopping time, we need to determine the pebble size, s_p . We consider the internal structure of pebble such that ice is coated on top of the silicates (“single-seed model” in the terminology of Schoonenberg & Ormel 2017). When pebbles cross the snowline and ice sublimates, it is assumed that the bare silicate core remains intact (Spadaccia et al. 2022; Houge et al. 2024, and see Saito & Sirono 2011; Hyodo et al. 2019 for case where pebbles disaggregate). Since the mass of the silicate core (mcore) in each pebble is unaffected by ice sublimation, the pebble number density can be obtained as $$n_{\text{p}}=\frac{{\rho }_{\text{p},\text{sil}}}{m_{\text{core}}},$$(13)

where ρ_p,sil is the volume density of silicate. From the number density, the particle mass is obtained, (ρ_p,sil + ρ_p,ice)/n_p, which, combined with its internal density, gives s_p : $$\frac{4\pi }{3}{s}_{\text{p}}^3=\frac{1}{{\rho }_{•,\text{p}}}\frac{{\rho }_{\text{p},\text{\hspace{0.17em}sil\hspace{0.17em}}}+{\rho }_{\text{p},\text{\hspace{0.17em}ice\hspace{0.17em}}}}{n_{\text{p}}}=\frac{m_{\text{core\hspace{0.17em}}}}{{\rho }_{•,\text{\hspace{0.17em}sil\hspace{0.17em}}}}\left(1+\frac{{\rho }_{\text{p},\text{\hspace{0.17em}ice\hspace{0.17em}}}}{{\rho }_{\text{p},\text{\hspace{0.17em}sil\hspace{0.17em}}}}\frac{{\rho }_{•,\text{\hspace{0.17em}sil\hspace{0.17em}}}}{{\rho }_{•,\text{\hspace{0.17em}ice\hspace{0.17em}}}}\right).$$(14)

In summary, the stopping time of the pebble is jointly determined by the densities of silicate and ice in the pebble compound. In the multifluid dust module, at each grid cell, the silicate and ice fluids share the same stopping time.

Apart from drift, pebbles also diffuse. Following Youdin & Lithwick (2007), the diffusivity of solids is related to the gas diffusivity through the Schmidt number Sc, $$D_{\text{p}}=D_{\text{g}}\times \text{Sc}=\frac{D_{\text{g}}}{1+{\tau }_{\text{s}}^2},$$(15)

where τ_s = Ωt_s is the dimensionless stopping time (Stokes number).

Throughout this work, we fix f_i/g = 0.4, which corresponds to a pebble-to-gas flux ratio of 0.8. Realistic disk evolution models invoking pebble growth suggest that f_i/g varies with e.g., disk metallicity, disk size and time (Birnstiel et al. 2010, 2012; Lambrechts & Johansen 2014; Ida et al. 2016; Schoonenberg et al. 2018). As shown in Schoonenberg et al. (2018), the pebble- to-gas flux ratio can readily reach unity at early times (<10⁵ yr) of the disk under solar metallicity. At later times, the flux ratio may drop a bit (Ida et al. 2016) but we do not expect this to qualitatively affect our results. The initial pebble size is set such that the Stokes number τ_s is τ_s,₀ at r₀. In all runs, τ_s,₀ = 0.03 is used, following Schoonenberg & Ormel (2017), which corresponds to a particle radius of s_p ≈ 2 cm in our disk model. With ice sublimation and vapor recondensation, the pebble size can decrease and increase following Equation (14). We do not include coagulation processes in our model.

2.4 Phase change

When pebbles drift across the snowline, the water ice component sublimates. This process transfers mass from ice to vapor, but also absorbs latent heat from the environment. Furthermore, the released water vapor, which has a higher mean molecular weight, will alter the pressure support thus the pressure scaleheight of the disk.

To properly capture the hydrodynamic and thermodynamic effect of pebble sublimation, we apply the phase change module developed in Wang et al. (2023). This module treats mass transfer and energy exchange during ice sublimation and vapor condensation self-consistently. The water vapor is treated as an additional tracer fluid, which advects with the gas flow and diffuses according to the diffusivity D_g .

Sublimation (condensation) occurs following rate expression (e.g., Ros & Johansen 2013; Schoonenberg & Ormel 2017). Within one timestep (Δt), the amount of vapor gained during sublimation (Δρ_vap)_subl or ice gained during condensation (Δρ_ice)_cond is $${\left(\mathrm{\Delta }{\rho }_{\text{vap\hspace{0.17em}}}\right)}_{\text{subl\hspace{0.17em}}}={\left(\mathrm{\Delta }{\rho }_{\text{p,ice\hspace{0.17em}}}\right)}_{\text{cond\hspace{0.17em}}}=\pi s_{\text{p}}^2{v}_{\text{th,vap\hspace{0.17em}}}{\rho }_{\text{vap\hspace{0.17em}}}{n}_{\text{p}}\left|1-\frac{P_{\text{eq\hspace{0.17em}}}}{P_{\text{vap\hspace{0.17em}}}}\right|\cdot \mathrm{\Delta }t,$$(16)

where n_p is the pebble number density, v_th,vap is the thermal velocity of vapor, P_vap is the vapor partial pressure, and P_eq is the saturation pressure of vapor as given by the Clausius–Clapeyron equation: $$P_{\text{eq}}=P_{\text{eq},0}\exp \left(-T_a/T\right).$$(17)

Here, P_eq,0 and T_a are constants specific to the species. For water, we take T_a = 6062 K and P_eq,₀ = 1.14 × 10¹³ g cm⁻¹s⁻² (Lichtenegger & Komle 1991). However, sublimation (condensation) stalls when the vapor equilibrium state given by Equation (17) is reached. Therefore, in the phase change module, the amount of material that is transferred between the phases is subject to the condition that vapor equilibrium is obtained whenever possible, namely, $$\begin{array}{l}\mathrm{\Delta }{\rho }_{\text{vap\hspace{0.17em}}}=\min \left[{\left(\mathrm{\Delta }{\rho }_{\text{vap\hspace{0.17em}}}\right)}_{\text{subl\hspace{0.17em}}},{\rho }_{\text{eq}}-{\rho }_{\text{vap}}\right],\\ \mathrm{\Delta }{\rho }_{\text{ice\hspace{0.17em}}}=\min \left[{\left(\mathrm{\Delta }{\rho }_{\text{ice\hspace{0.17em}}}\right)}_{\text{cond\hspace{0.17em}}},{\rho }_{\text{vap\hspace{0.17em}}}-{\rho }_{\text{eq}}\right],\end{array}$$(18)

where ρ_eq is the vapor density corresponding to the saturation pressure. In this way, at the disk midplane around the snowline, vapor tends to follow the saturation pressure, while at the disk atmosphere, vapor tends to be super-saturated due to the slow condensation rate (Equation (16)) and can be freely transported by diffusion and advection.

Latent heat is defined as the enthalpy difference between the vapor and the ice. Supposing that δρ is the mass transfer from ice to vapor during sublimation, an additional energy δE = L_heatδρ should be subtracted from the total internal energy, while the same amount of energy will be released if condensation occurs. For water ice, we take L_heat = 2.75 × 10¹⁰ erg g⁻¹ (Fray & Schmitt 2009). For the description of the phase change module and its detailed implementation in Athena++, we refer to Wang et al. (2023).

The ideal equation of state of gas is adopted in this work, $$P_{\text{g}}={\rho }_{\text{g}}\frac{k_{\text{B}}T}{\mu m_{\text{p}}}.$$(19)

Once vapor is released, we change the mean molecular weight of the gas accordingly following, $$\frac{1}{\mu }=\frac{f_{\text{v}}}{{\mu }_{{\text{H}}_2\text{O}}}+\frac{1-f_{\text{v}}}{{\mu }_{\text{xy}}},$$(20)

where f_v is the vapor mass fraction, and µ_H2O = 18 and µ_xy = 2.34 are the mean molecular weight of water and H-He mixture, respectively.

2.5 Two-stream radiation transfer

In this work, we focus on the steady state solution of the snowline region rather than dynamical and time-dependent behaviors. Thus, we opt to determine the temperature structure using an equilibrium solution obtained with the dissipation profile. Following Mori et al. (2019), we define the total dissipation profile as $$q_{ɀ}=q_{\text{irr}}+q_{\text{vis}}+q_{\text{latent\hspace{0.17em}}}+q_{\text{diff}},$$(21)

where q_irr, q_vis, q_latent, and q_diff are the heating terms attributed to stellar irradiation, viscous dissipation, latent heat exchange, and heat diffusion, respectively. Here, we implicitly adopt the two-stream assumption, where the incoming stellar irradiation and disk thermal emission are well separated in the visible and infrared band, respectively, so that we can combine the corresponding heating sources independently.

Following Calvet et al. (1991), the heating rate contributed by stellar irradiation, q_irr , is $$q_{\text{irr}}=E_0\rho {\kappa }_{\text{vi}}\left[\exp \left(-\frac{{\tau }_{\text{vi}(ɀ)}}{{\mu }_0}\right)+\exp \left(-\frac{{\tau }_{\text{vi}}(-\infty )-{\tau }_{\text{vi}}(ɀ)}{{\mu }_0}\right)\right].$$(22)

Here, E₀ = L_★/(8πr²) is the stellar irradiation flux and κ_vi, τ_vi are the opacity and optical depth for visible light, respectively: $${\tau }_{\text{vi}}={\displaystyle {\int }_{ɀ}^{+\infty }\rho }{\kappa }_{\text{vi}}\text{\hspace{0.17em}d}{ɀ}^{\prime }.$$(23)

Also, µ0 is the cosine of the incident angle of the stellar irradiation, $${\mu }_0=r\frac{\text{\hspace{0.17em}d}}{\text{\hspace{0.17em}d}r}\left(\frac{H_{\text{p}}}{r}\right),$$(24)

and H_p is the height of the photosphere of the disk (Chiang & Goldreich 1997). We took H_p = 4H_g in our modeling. In principle, ray tracing is needed to calculate the irradiative heating accurately. However, this simplification will not affect the overall results since the ice-sublimating regions are always well within optically thick limit for stellar irradiation in all runs.

For active disks, where viscous dissipation is deposited in the midplane, the viscous heating term q_vis (e.g., Armitage 2020) follows $$q_{\text{vis}}=\frac{9}{4}\rho \nu {\mathrm{\Omega }}^2.$$(25)

To include the thermal effect of ice sublimation (condensation), we model the latent heat absorption (release) as an additional heating term defined as $$q_{\text{latent\hspace{0.17em}}}=-\frac{\text{d}{\rho }_{\text{vap\hspace{0.17em}}}}{\text{d}t}\cdot L_{\text{latent\hspace{0.17em}}}.$$(26)

The change of vapor mass dρ_vap/dt is measured in the simulation, as determined by the phase change module.

Lastly, in optically thick regions, unlike optically thin regions where vertical transfer is usually quicker than the radial component due to the small aspect ratio of the disk, heat diffusion works both vertically and radially. We aim to complete our radiative transfer approach by including the radial heat diffusion only in the optically thick region. Therefore, following Xu & Kunz (2021), $$q_{\text{diff\hspace{0.17em}}}=\nabla \cdot \left[{\kappa }_{\text{heat\hspace{0.17em}}}\nabla T\exp \left(-1/{\tau }_{\text{R}}\right)\right],$$(27)

and, $${\kappa }_{\text{heat\hspace{0.17em}}}=\frac{16{\sigma }_{\text{SB}}{T}^3}{3{\kappa }_{\text{R}}\rho }.$$(28)

Here κ_heat is the heat conduction coefficient evaluating the efficiency of energy transport by radiative diffusion (e.g., Kippenhahn et al. 2013), κ_R the Rosseland mean opacity and τ_R the corresponding optical depth. The exponentially diminishing term exp (−1/τ_R) suppresses the heat diffusion in optically thin regions.

Given the total dissipation profile, q_z , the equilibrium solution of the temperature can be expressed as $$T_{\text{eq}}(ɀ)=T_{\text{eff}}{\left(\frac{3}{4}{\tau }_{\text{eff}}+\frac{\sqrt{3}}{4}+\frac{q_{ɀ}}{4\rho {\kappa }_R{ℱ}_{+\infty }}\right)}^{1/4}.$$(29)

We refer the detailed derivation of Equation (29) to previous works in the literature (e.g., Hubeny 1990; Mori et al. 2019), but we briefly describe the meaning of the expression here. First, ℱ_+∞ is the radiative flux leaving the upper surface and T_eff is the corresponding effective temperature: $${ℱ}_{+\infty }={\displaystyle {\int }_0^{+\infty }{q}_{ɀ}}\text{d}ɀ;T_{\text{eff\hspace{0.17em}}}={\left(F_{+\infty }/{\sigma }_{\text{SB}}\right)}^{1/4}.$$(30)

Then, τ_eff is the flux-weighted effective optical depth, $${\tau }_{\text{eff}}(ɀ)=\frac{1}{{ℱ}_{+\infty }}{\displaystyle {\int }_{ɀ}^{+\infty }\rho }{\kappa }_{\text{R}}ℱ\left({ɀ}^{\prime }\right)d{ɀ}^{\prime },$$(31)

where $ℱ(z)={\displaystyle {\int }_0^z{q}_z}dz$. If a high amount of heating is deposited in optically thick regions, τ_eff will dominated Equation (29). In contrast, in the optically thin regions, the temperature is determined by a local balance of heating and cooling terms (Eq. (29), third term) and the radiation field sets the background equilibrium temperature (Eq. (29), second term).

In the simulation, in order to stabilize the system, the temperature is gradually relaxed by introducing a thermal relaxation parameter β: $$\frac{\text{d}T}{\text{\hspace{0.17em}d}t}=\frac{T_{\text{eq}}-T_0}{\beta {\mathrm{\Omega }}^{-1}},$$(32)

where T₀ is the temperature before the thermal relaxation. Although this relaxation breaks the equilibrium solution obtained with the two-stream radiation transfer method, the steady state solution will remain consistent, since dT /dt ≈ 0. We choose β = 50 in our simulation to prevent numerical oscillations.

2.6 Simulation cases

Our aim is to investigate the impact of qualitatively distinct disk thermal structures on the characteristics of the snowline region. Consequently, we define three disks:

• Passive disks. The temperature is determined by how deep the stellar irradiation (q_irr ) penetrates while the contribution from viscous heating is negligible. The disk is optically thick for cooling. Passive disks feature a hot upper layer and a cold midplane.

• Active Disks. Here the viscous heating (q_vis) dominates, which, combined with the optically thick environment that prevents cooling, renders the temperature higher near the midplane.

• Isothermal Disks. The main heating source is also the stellar irradiation as the passive disks. However, the disk is optically thin for cooling.

These three cases represent the main snowline thermal morphology investigated in this work.

Because we aim to directly compare the simulations to previous 1D calculations (Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017; Hyodo et al. 2021), we set up the disk model by fixing α = 3 × 10⁻³ and Ṁ_acc = 10⁻⁸ M_⊙ yr⁻¹, following Schoonenberg & Ormel (2017). In addition, we anchor the midplane snowline location (r_snow,_mid) at ≈2.1 au before the injection of icy pebbles (which will shift r_snow when sublimating) by tuning the stellar luminosity and opacity (see below). The location of the snowline is approximated by equating the vapor saturation pressure to the expected vapor pressure if all incoming ice is sublimated (Schoonenberg & Ormel 2017). $$P_{\text{eq},0}\exp \left[-\frac{T_a}{T}\right]=\frac{k_{B}T}{\mu }{\rho }_{\text{g},\mathrm{mid}}{f}_{\text{i}/\text{g}},$$(33)

where ρ_g,mid is the midplane gas density. In practice, given a disk viscosity profile and mass flux, the surface density profile can be obtained. Fixing the surface density at a certain radius, the gravity-balanced vertical density and temperature distribution can be iteratively solved. This serves as the initial condition of the simulation and supplies the initial estimate of r_snow .

Under this simplified disk model, the main parameters that determine the initial thermal structure are κ_vi , κ_R, and L_⋆ , which are all highly uncertain among different disk environments (e.g., Birnstiel et al. 2018). To simplify, we use constant Rosseland- mean opacity with respect to the density of non-condensable gas (H and He) τ_R = ∫κ_Rρ_H,Hedz (Table 1). Also, we use the same visual band opacity κ_vi = 10 cm²g⁻¹ among all runs. This value corresponds approximately to the opacity of a grain size distribution with the maximum grain size (a_max) around sub-mm and a power-law of −3.5 (Mathis et al. 1977), as calculated from the DSHARP opacity model (Birnstiel et al. 2018).

By varying κ_R and L_⋆, we were able to construct disks with qualitatively-distinct thermal morphology and their initial midplane snowlines at 2.1 au. As shown in Fig. 1, the permitted parameter space can be divided into three regimes: “passive disk” when L_⋆ is high, “active disk” when κ_R is large and “vertically-isothermal” disk when κ_R is small. The chosen parameters of simulation runs are denoted by black dots in Fig. 1 and summarized in Table 1.

As test cases to validate our model, we also conduct 1D snowline simulations (Table 1). The 1D/SO17 run is comparable to the “simple, single-seed” model in Schoonenberg & Ormel (2017) except that the increase of mean molecular weight with vapor injection is included here. The 1D-hydro run has the gas motion solved hydrodynamically, yielding distinct results, which is discussed in Appendix B and Sect. 6.1.

Fig. 1 Initial midplane snowline location resulting from the two-stream radiation transport model (Equation (29)) as function of Rosseland mean opacity κR and stellar luminosity L⋆ . Here the optical opacity is fixed to κvi = 10 cm2g−1. The permitted parameter space is divided into three regimes (white dashed lines, for illustration purpose only): “passive disk” when L⋆ is high, “active disk” when κR is high and “vertically-isothermal” disk when κR is small. The dotted lines are temperature contours of the midplane at 2.1 au, ranging from 100 K to 200 K in 10 K increments. The black solid line represents the parameter combinations that lead to rsnow,mid = 2.1 au and the black dots denote the parameter combinations used in this work, as summarized in Table 1. Hydrodynamic and thermodynamic effects will move the snowline location in the simulation.

3 Results

In this section, we present the results of the 2D disk snowline models as listed in Table 1 and introduced in Sect. 2.6. We first present the overall simulation setup and demonstrate that the simulations reached a steady state based on the measured mass fluxes (Sect. 3.1). Then we introduce the temperature structure, density structure and water delivery patterns (Sect. 3.2). In Sect. 3.3, we illustrate the effects of different snowline morphologies (active, passive, isothermal) have on the pile-up of solids at the snowline. We identified a “water cycle” that enhances the pile-up of water in the snowline region in active disks. The effects of latent heat will be discussed in Sect. 4.

3.1 Simulation setup and verification

The midplane snowline is initially anchored at 2.1 au (Sect. 2.6). We set up the simulation in a spherical-polar (R-θ) coordinate system with the star centered at R = 0. The simulation domain covers R ∈ [1.0,4.0] au and θ ∈ [1.3,π/2]. We describe our models (Sect. 2) and conduct analysis of the data in a cylindrical coordinate (r-z) for convenience. At the reference position r₀ = 3.0 au the simulation domain in θ corresponds to ≈7H_g. The grid resolution is N_r × N_θ = 320 × 64. To improve the resolution near the midplane region, where most of the material resides, we adopt a uniform spacing for θ^1/3 (e.g., Ormel et al. 2015). This increases the grid resolution near the midplane by a factor of ≈3.

Simulations are initialized with non-condensible gas (H-He) and relaxed to a steady flow pattern. Then pebbles are released from the outer boundary and sublimate at the iceline region to enrich the inner disk with vapor. We monitor the ice and vapor fluxes (ℱ_ice and ℱ_vap) and evolve the system until a steady state is reached. For all simulations, a steady state is reached after about 2 × 10⁵ Ω⁻¹. To illustrate the realization of the steady state, we plot the vertically-integrated radial mass flux of model passive in Fig. 4c. Exterior to the iceline the ice flux ℱ_ice equals the imposed f_i/gṀ_acc, until sublimation starts at around 2.5 au where the vapor flux becomes non-zero. In the snowline region, where phase change processes operate, the total water flux ℱ_water = ℱ_ice + ℱ_vap is conserved and equals the imposed ice flux. The silicate (ℱ_sil) and H-He gas flux (ℱ_xy) are also plotted in Fig. 4c and shown to be conserved across the simulation domain.

3.2 2D steady state structure and characteristic outputs of the solid pile-up

To illustrate the qualitative picture of the 2D snowline, we present the radial(r)-vertical(z) slices of the steady state density distribution in Fig. 2 and the corresponding temperature structure in Fig. 3. In these figures the green lines are drawn as the boundary ρ_p,ice/ρ_g = 10⁻², which represent the surface where ice is exhausted in the disk (we call it “snowline” here- after¹). In Fig. 2, the region enclosed by the green line is filled with blue contours that represent the ice density (normalized to the midplane gas density at r₀ of the unperturbed disk, i.e., ρ₀ ≈ 3.4 × 10⁻¹¹ cm g⁻³), while the rest is filled with red contours that represent the vapor density. To further visualize the water transport process, we plot ice (blue) and vapor (pink) mass fluxes density (in the unit of ρv) with streamlines, the thickness of which denote the strength of the mass flux. The mass flux density includes both the advective and diffusive components.

Consistent with previous 1D results (e.g., Schoonenberg & Ormel 2017), Fig. 2 demonstrates the characteristic ice density bump exterior to the snowline in all runs. This bump is more clearly illustrated in Fig. 4a, where the surface density profile of the passive run is presented. The features observed in this run are representative to all simulations runs. As pebbles drift in, the initially-equal surface density of ice and silicate (ζ = 0.5) gradually separates due to the pile-up of ice. The pile-up in solids just exterior to the snowline is mostly contributed by ice rather than silicate (Fig. 4a). Eventually all ice sublimates into vapor, which then follows the accretion flow of the gas. Therefore the surface density of vapor matches the expected viscous solution of gas in steady state (grey line), while silicate pebbles keep drifting towards the star at high velocity (typically ten times higher than the gas accretion). The corresponding midplane solid-to-gas ratio – an indicator of planetesimal formation (e.g., Johansen & Youdin 2007) – is plotted in Fig. 4b. To quantify the properties of the solid pile-up and compare its strength among different simulation runs, we defined the following indicators:

1. The location of the midplane snowline r_snow,mid as defined in Equation (33). While r_snow,mid is initially anchored to 2.1 au, it could change after ice injection.

2. The peak midplane solid-to-gas ratio ξ_pk = max [(ρ_p,ice + ρ_p,sil)/ρ_g] and its corresponding location r_pk.

3. The full-width-half-maximum (FWHM) of the solid pile-up as measured from the dust-to-gas ratio profile.

4. The ice mass (M_ice) stored in the solid pile-up. We integrate the ice surface density from the point where Σ_ice/Σ_sil > 1.5 inwards (blue-shaded region in Fig. 4a).

The positions of r_pk, r_snow,mid, and the FWHM of the solid pileup are illustrated in Fig. 4b and the value of these characteristic indicators are summarized in Table 1 for all runs.

In previous 1D studies (Ros & Johansen 2013; Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017; Hyodo et al. 2021), the solid pile-up is understood as a result of outward diffusion of vapor and recondensation. To examine this picture, the vertically-integrated radial mass flux of passive, active and active-nqL are plotted in Fig. 5 (upper panel). The vapor flux is further split into diffusive and advective components. Consistent with the flux density of vapor just exterior to the snowline shown in Fig. 2, we observe a strong outward flux of vapor in all three cases, responsible for the pile-up. In the passive run, the total outward vapor flux is dominated by diffusion. However, in the active disk runs (active-nqL and active), the contribution of the advective flux is comparable to that of the diffusive flux. This implies that there must be a strong outflow near the snowline region, which is confirmed by the midplane radial gas velocity profile in Fig. 5 (lower panel). While this radial outflow is automatically captured by the hydrodynamic simulation, 1D models relying on integrated versions of the transport equations (e.g., Schoonenberg & Ormel 2017) ignore it and therefore underestimate the pile-up strength (see the discussion in Appendix B). In Sect. 5.1, we present an analysis of this radial outflow to explain its origin.

Aside from the snowline region, water transport in general is characterized by the inward drift of ice from the outer boundary and the accretion of water vapor at the inner boundary (Fig. 2). Focusing on the disk region inside the snowline, where the mass flux is purely contributed by advection (see Fig. 5, upper panel), the streamlines of vapor flux density (pink) reveal outflow at the midplane and accretion at the upper regions (≥1 H_g). This characteristic flow pattern is consistent with previous studies of 2D viscous disk with pure H-He gas (Takeuchi & Lin 2002; Ciesla 2009). The preferential gas accretion in the upper layer is due to the steep radial density gradient at the midplane, resulting in a net gain of angular momentum due to the greater viscous torque exerted by the faster-rotating inner disk compared to the slower-rotating outer disk (e.g., Takeuchi & Lin 2002). Combined with the outflow in the midplane and the outward diffusion of vapor near the pile-up region, the simulations demonstrate that vapor leaves the snowline region and is accreted easier from higher altitude (see also Sect. 5.2).

Fig. 2 Steady state density structure of all simulated disk (as listed in Table 1). The dark green line denotes the boundary where ρp,ice/ρg = 10−2. The region enclosed by the green line filled with blue shading represents the ice density normalized to ρ0 ≈ 3.4 × 10−11 cm g−3 (the midplane gas density at r0 of the unperturbed disk), while the region interior to the snowline indicates the vapor density with red shading. The ice (blue) and vapor (pink) mass fluxes are plotted with streamlines, whose thickness denotes the magnitude of the flux normalized to ρ0cs,0. The mass flows in the vertical direction, which are especially vigorous in the active-nqL and active runs, contribute towards an increased pile-up of ice in the snowline region.

3.3 Effect of snowline morphology

We go on to analyze the effect of the snowline morphology, first focusing on the simulation runs neglecting latent heat effects (active-nqL, passive-nqL and iso-nqL). The temperature structures (Fig. 3) of the passive and isothermal disks looks similar, resulting in a similar shape of the snowline, which shifts inwards more at the midplane due to the higher ice density. In contrast, the active-nqL run harbors a hotter midplane, resulting in an inward protrusion of the snowline above the midplane, similar to the findings by Oka et al. (2011).

Different snowline morphologies lead to distinct pileup properties. From Table 1, the active-nqL run stores a significantly higher ice mass (M_ice) in the pile-up region. The active-nqL run also has a narrower pile-up (FWHM), which results in a higher peak dust-to-gas ratio, ξ_pk. The smaller FWHM can be understood since the temperature gradient in the active-nqL run is larger across the snowline, as reflected by the closer-spaced contours in Fig. 3. As a result, the ice sublimation takes place in a narrower region.

Furthermore, from the ice flux density streamlines in Fig. 2a, we see that the active-nqL run is characterized by strong pebble settling in the pile-up region, a feature that is absent in the passive and isothermal disks. In all runs, ice sublimates mainly at the midplane near the snowline. The released vapor then spread upwards (mostly by diffusion). Since it is cooler in the upper layers in the active disk, vapor recondenses on pebbles, which then settle back to the midplane. This process is less efficient in the passive and isothermal disks since the temperature is more uniform in the vertical direction. As pebbles settle, they sublimate again, which leads to a “cycle” of pebble growth and settling. In a steady state, we observe that the Stokes number of pebbles increases to ≈0.1–0.2 near the upper boundary of the snowline, in contrast to the nominal value of ≈0.03 at the midplane, which leads to a strong downward flux of ice, as seen in Fig. 2a.

As discussed in Sect. 3.2, vapor escapes the snowline region predominantly from higher altitude. As a result, the recondensation of vapor in the upper layers of active disks effectively confines water to the snowline region. On the other hand, in the passive and isothermal disks, vapor can freely flow to the inner disk after it has spread to the upper layers. This water cycle observed in the active disk results in higher amounts of ice stored in the snowline region (Table 1). A further analysis will be presented in Sect. 5.2.

Fig. 3 Steady state temperature structure of all the simulation runs (Table 1). The color shows the temperature structure with contours highlighting the range between 150 and 200 K in increments of 5 K. The dashed lines show the Rosseland mean optical depth τR from 0.1 to 10.0, displaying the optically-thick or optically thin nature of the cooling. In panel (e) and (f), the dashed lines are absent since the disk is very optically thin.

3.4 Effects on the rotational velocity

With vapor injection, the pressure support near the snowline will be altered due to the sudden increase of density and mean molecular weight, therefore changing the azimuthal velocity (v_ϕ). The strength of the radial pressure gradient is usually expressed through (Nakagawa et al. 1986): $$\eta \equiv -\frac{1}{2}\frac{c_{\text{s}}^2}{v_{\text{k}}^2}\frac{\partial \log P_{\text{g}}}{\partial \log r}\simeq \frac{v_{\text{hw}}}{v_{\text{k}}},$$(34)

where v_hw = v_k − v_ϕ ≃ ηv_k denotes the deviation of the gas motion from Keplerian motion, i.e., the headwind (e.g., Armitage 2020). We measure η at the midplane of the disk before vapor injection and after the steady state is reached (Fig. 6). Compared to disks without vapor injection (“unperturbed”), η decreases for all runs inside the snowline. In addition, η exhibits a dip right at the snowline location, r_snow,mid . To understand this feature, we expand the pressure gradient in Equation (34) and write η as $$\eta =-\frac{1}{2v_{\text{k}}^2}\frac{k_{\text{B}}T}{\mu m_{\text{p}}}\frac{\text{\hspace{0.17em}d}}{\text{\hspace{0.17em}d}\log r}\left(\log {\rho }_{\text{g}}+\log T-\log \mu \right).$$(35)

The negative sign of the µ-gradient term arises due to the decreasing pressure support from a higher mean molecular weight gas.

Inside the snowline, the decrease in the headwind is fully caused by the overall increase of µ, i.e., the prefactor in Equation (35). Given that the water vapor fraction is f_v = f_i/g/(1 + f_i/g) = 0.4/1.4, µ will increase from µ_xy = 2.34 to µ = 3.11 (Equation (20)), consistent with the shifts of η in Fig. 6. At the snowline, the dip in η results from the steep gradient in the mean molecular weight (−d log µ/d log r), which decreases pressure support, although the vapor injection (d log ρ_g /d log r) from ice sublimation counteracts this effect. The dip is most pronounced in the active-nqL run where the vapor fraction increases steepest in the snowline region. Surprisingly, the active run exhibits a similar level of decrease in η as the active-nqL run, despite its smaller gradient in µ. The additional decrease in η in the active run originates from the significant flattening of the temperature profile (log T ) by latent heat cooling, which will be discussed in Sect. 4.

In summary, compared to pure H-He disk without vapor injection, the headwind generally decreases by a factor of 0.6– 0.7 at the snowline, which is significant enough to boost pebble accretion rates and promote SI (see Sect. 6.2).

Fig. 4 Radial profiles of the passive run. (a): vertically-integrated surface density of total gas, ice, silicate and vapor. The grey line shows the expected vapor surface density from the unperturbed viscous solution (Equation (2)), where fi/g is the ice-to-gas flux ratio. The blue shade denotes the region that contains Mice. (b): the solid/vapor-to-gas ratio in the midplane. The positions of peak solid-to-gas ratio (rpk) and the midplane snowline (rsnow,mid) are highlighted; (c): vertically- integrated radial mass fluxes of ice, vapor, water, H-He and silicate, where ℱwater = ℱvap + ℱice. The dashed lines denote the imposed gas (Ṁacc) and ice ( fi/g Ṁacc) fluxes, respectively.

4 Effect of latent heat exchange

Latent heat exchange associated with ice sublimation and vapor recondensation has not been investigated explicitly in disks with hydro simulations before (see the discussion of Owen 2020). To examine its impact on the thermal structure of the disk and the resulting solid pile-up, we run a series of simulations of different snowline morphology incorporating q_latent, the latent heat exchange term. All other model parameters are kept the same (see Table 1) to enable a direct comparison.

4.1 Active disks

In the active disk, the temperature structure near the snowline is significantly altered. Unlike the smooth, gradual temperature transition across the snowline seen in active-nqL (see the evenly-spaced contours in Fig. 3a), the active run (Fig. 3b) reveals a slowly-varying temperature “plateau” in the ice sublimation region, followed by a sharp temperature increase across the snowline (marked by green line). This is further illustrated in Fig. 7b, where a jump in the midplane temperature (solid line) is clearly noticeable. The plateau arises from the latent heat absorption during ice sublimation. A maximum temperature difference of ≈40 K is observed when q_latent is included, highlighting the significant cooling effect due to latent heat absorption.

With the thermal structure altered by ice sublimation, the properties of the pile-up are also affected. First, the disk snowline is pushed inwards (see the green lines in Figs. 3a, b), resulting in a shift in the location of the pile-up (Fig. 7a). Second, both the width and the amplitude of the pile-up are modified. As shown in Fig. 7c, the temperature plateau created by latent heat cooling from ice sublimation broadens the snowline region. Consequently, the pile-up becomes around two-times more extended, while the peak solid-to-gas ratio is reduced by a factor of ≈1.6 (Table 1). The extended snowline region in theactiverun is also evident from Fig. 2b, where pebble settling happens over a wider region compared to the active-nqL run.

Comparing the temperature contours in Figs. 3a and b, we find that the latent heat exchange flattens the temperature gradient both in the vertical and radial directions. In Fig. 7b, we present the temperature profiles at various altitudes by selecting slices in θ-direction. For the midplane and ɀ(r_snow) = 0.04 au slices, the temperature is nearly uniform in the snowline region of theactiverun, whereas a distinct vertical temperature gradient is present in the active-nqL run. The more isothermal vertical temperature structure results from the combination of latent heat absorption due to ice sublimation at the hotter midplane and latent heat release from vapor recondensation at the cooler upper layers. As a result, ice settling is mitigated in the active run compared to the active-nqL run (Fig. 2, blue streamlines). However, accounting for latent heat exchange does not overturn the fundamental thermal stratification, i.e., the midplane remains hotter than the upper layers (Fig. 3b). As a result, the active disk exhibits a water cycle similar to that in the active-nqL disk (Fig. 2) and the total ice mass (M_ice) stored in the snowline region is not significantly influenced by latent heat exchange (Table 1).

Fig. 5 Upper panels: vertically-integrated radial mass flux of vapor. For three runs (passive,activeand active-nqL), the diffusive, advective and total vapor flux are plotted. The dotted lines denote zero-flux levels. The blue arrows indicate the position of the midplane snowline (rsnow,mid). The advective flux contributes significantly to the total vapor flux in theactivedisk runs (active and active-nqL), while it is negligible in the passive run. Lower panels: midplane gas radial velocity and the contributions from different terms following Equation (36). The thick grey line represents the measured value from the simulation. The T2 term signifies the dominant contribution of the radial viscous stress to the radial velocity, while the contributions from radial advection (T1 ) and ɀ-direction (Tɀ) are negligible.

Fig. 6 Azimuthal velocity deviation at the midplane , expressed in terms of the radial pressure gradient parameter η (Equation (35)). The black dashed and solid lines show η before vapor injection and after the steady state is reached, respectively. The red lines show the vapor fraction corresponding to the y-axis on the right. The blue arrows indicate the position of the midplane snowline (rsnow,mid). A reduction of η at the snowline region strengthens solid pileups and is conducive to planetesimal formation by streaming instability and planet formation by pebble accretion.

Fig. 7 Comparisons of radial profiles of the active disks. For all three panels, the solid line represents the active disk while the dashed line represents active-nqL disk. (a) The surface density profiles. As in Fig. 4, the expected vapor surface density fi/g Ṁacc /(3πν) is plotted in grey. (b) The temperature profile of slices in θ-direction at the midplane and higher altitudes. The slices above the midplane pass through a certain height ɀ at the location of snowline (rsnow) of the active run. (c) The solid-to-gas ratio (black) and vapor-to-gas ratio (red) at the midplane in steady state.

4.2 Passive and isothermal disks

Although q_latent reshapes the snowline structures in active disks, its impact on passive and isothermal disks is significantly less pronounced, as evidenced by the characteristic indicators of the pile-up (Table 1). This is clearly seen by comparing the steady state thermal structure of runs with and without q_latent (Fig. 3). In active disks, the temperature contours are radially stretched and compressed over the range of around 1.7−2.3 au at all altitudes within the phase change region. However, the temperature structure in passive is nearly identical to that in passive-nqL . In the isothermal disks, the temperature contours are only slightly bent near the snowline, implying a much weaker and more localized effect of latent heat.

These different effects of latent heat on the thermal structure can be explained by their distinct leading terms contributing to the temperature from Equation (29). For passive disks in our simulations, due to the large L_⋆ chosen (Table 1), the second (constant) term dominates, representing the radiationequilibrium temperature determined by the stellar radiation field. Consequently, the contribution from q_latent is negligible. Similarly, in isothermal disks, since it is very optically thin (κ_R very low), the third term in Equation (29), which represents the immediate balance between heating and cooling due to q_ɀ , matters. This energy balance is local, therefore, the latent heat effects are also localized.

In contrast, in active disks, where the snowline region is optically-thick (see the dotted lines in Fig. 3) and L_⋆ is low, the “blanketing” effect, related to the effective, namely, flux- weighted, optical depth τ_eff (first term in Equation (29)), is dominant. Energy deposited deeper within the optically-thick disk heats more efficiently, as it takes longer for the energy to escape (e.g., Chandrasekhar 1935; Mori et al. 2021). Since ice mainly sublimates at the midplane, latent heat absorption directly counteracts viscous heating, reducing the net heating (q_ɀ ≪ q_vis), which significantly mitigates the “blanketing” effect.

To summarize, the effect of latent heat is most pronounced in disks where the blanketing effect is strong while it diminishes when cooling is efficient or when the disk is in radiation equilibrium.

5 Analysis

5.1 Gas dynamics near the snowline

As discussed in Sect. 3.2, there is a strong outflow of gas near the snowline. This is clearly shown in Fig. 5, where the midplane gas velocity near the snowline is larger than other unperturbed region (e.g., ≈10 times larger in the active-nqL run). To understand it, we derive the expression for radial gas velocity in steady state, starting from the continuity and Navier- Stokes equations (considering angular momentum conservation; e.g., Balbus & Papaloizou 1999; Jacquet 2013). In cylindrical coordinate (r-ϕ-ɀ), the radial velocity of gas is given by $$\begin{array}{c}{\nu }_r=-\frac{2}{r\rho {\nu }_{\text{k}}}\frac{\partial }{\partial r}[\underset{T_1}{\underbrace{r^2\rho v_r\delta v_{\varphi }}}\underset{T_2}{\underbrace{+\frac{3}{2}r\rho v_{\text{k}}v-r^3\rho v\frac{\partial }{\partial r}\left(\frac{\delta v_{\varphi }}{r}\right)}}]\\ \underset{T_{ɀ}}{\underbrace{-\frac{2}{\rho v_{\text{k}}}\frac{\partial }{\partial ɀ}\left[r\rho \delta v_{\varphi }{v}_{ɀ}-r\rho v\frac{\partial \delta v_{\varphi }}{\partial ɀ}\right]}}\end{array},$$(36)

where v_k represents the Keplerian velocity, δv_ϕ = v_ϕ − v_k is the deviation of azimuthal velocity from Keplerian and v_ɀ is the vertical velocity component. The steady-state radial velocity v_r results from the angular momentum (AM) transport by advection and viscous stress in both radial and vertical directions. The terms T₁ and T₂ represent the contributions from the radial advection and viscous stress, respectively, while T_ɀ encapsulates the AM transport from the ɀ-direction as a whole.

Evaluating Equation (36), Fig. 5 shows the gas radial velocity v_r,mid in the midplane, where the outflow is the strongest. The contributions from each term are plotted and the sum of these terms matches the measured radial velocity, validating the decomposition. Across the snowline, there is a strong outward motion (positive) of gas followed by a sharp transition to inflow (negative). The term T₂ dominates over other terms in v_r,mid, indicating that both outflow and inflow are primarily driven by viscous stress. This process can be understood as follows: as pebbles drift across the snowline and sublimate, the released vapor viscously spreads both outward and inward. In the absence of vapor injection, midplane gas motion is already characterized by outflow driven by the density gradient (see Sect. 3.2). The injection of vapor into the system further steepens the initial radial density profile, exacerbating the outflow. Compared with the background outflow inside the snowline, as identified in previous literature (Takeuchi & Lin 2002; Ciesla 2009), vapor injection induces >10 times stronger outflow at the snowline (Fig. 5, lower panel). Furthermore, the strength of the viscous spreading is determined by the density of the deposited vapor. In the active disk, the latent heat cooling creates a temperature plateau (Fig. 7b), leading to a more gradual vapor deposition near the snowline (Fig. 7a), and consequently a weaker outflow compared to the active-nqL disk. Similarly, the inflow is stronger in the active disk because there is a steeper increase of temperature interior to r_snow,mid, where ice is fully consumed and the latent heat cooling abruptly ceases. In summary, radial outflow from viscous spreading of vapor is commonly seen near the snowline.

5.2 Water-cycle enhances solid pile-up

In Sect. 3.3, we identify a water-cycle that traps a large amount of ice in the snowline region in active disks. In other words, water parcels take longer time to escape the snowline region in active disks. To examine this finding, we perform Lagrangian trajectory integrations of water particles in the active-nqL and passive-nqL runs. The trajectories of water particles – both in ice form as well as vapor – is modelled with a Monte Carlo method, which includes a stochastic component, representing particle diffusion, as well as a systematic component, representing the steady-state advective velocity fields (Ciesla 2010; Krijt et al. 2016). The simulation domain is mirrored to allow particles to travel into the lower hemisphere of the disk. At any point, the phase of the water parcel is determined from the sublimation and condensation rates, which represent the probability to either transfer from ice to vapor or the other way around. See Appendix A for a detailed description.

For each disk, 1000 particles (N_p = 1000) are released from the outer boundary and integrated for 10⁵ Ω⁻¹ or until the time where they get accreted (r < 1.1 au). The diffusion timescale for a water particle to spread over the gas scale height is $t_{\text{dif}}\sim H_{\text{g}}^2/D_{\text{g}}\approx 300{\mathrm{\Omega }}^{-1}$. Therefore the integration time is sufficient for the particles to reach a diffusion-advection equilibrium state in the vertical direction. And as the radial transport time much exceeds the vertical transport time, trajectories are characterized by a large amount of vertical movement. After obtaining the trajectories, we count the time that particles spent in each grid cell (t_spend) and plot its distribution in Figs. 8a,b. We observe a significant increase of t_spend near the snowline region, especially concentrated around the peak position of the pile-up (r_pk). This consistency implies that the pile-up is indeed driven by longer escape time of particles.

Clearly, water particles spend significantly more time in the snowline region of the active-nqL disk compared to the passive-nqL disk, explaining the stronger pile-up of M_ice in active disks. In other words, water is trapped in the active disks. Efficient transport of water away from the snowline region primarily occurs when water diffuses to upper layers (ɀ ≈ 0.1 au), where it is carried by the advective accretion flow (shown by the long red arrows in Figs. 8a,b). Near the midplane, however, only “lucky” particles can escape the snowline region by radial diffusion. This inefficiency arises because (1) vapor motion near the midplane (ɀ ≈ 0–0.05 au) is predominantly governed by diffusion, as depicted by the tiny red arrows in Figs. 8a,b; and (2) strong outflows along the snowline surface (thick grey lines), especially in the active-nqL run, create a significant barrier to the inward motion of water particles.

Moreover, due to the cold upper layer, strong ice settling (blue arrows in Figs. 8a,b), as part of the water-cycle, confines water much closer to the midplane in the active-nqL run. This is reflected in the root-mean-square of the particles’ vertical positions $\left(\sqrt{\langle {ɀ}^2\rangle }\right)$ (indicated by black dashed lines in Figs. 8a,b), which represents the scaleheight of water parcels. In Fig. 8a, it is also evident that at the radial location where t_spend is particularly large at the midplane, the water scaleheight tends to be compressed. As a result, the key upper-layer channel to escape the snowline region is suppressed in the active-nqL run, leading to a stronger pile-up of ice. The suppression of the upper-layer escape channel also explains why an outward, vertically integrated advective flux of vapor, comparable in magnitude to the diffusive flux, is observed in the active disk runs (active and active-nqL) but not in the passive run (Fig. 5 and Sect. 3.2). While all runs exhibit outflows that drive the outward advective flux of vapor at the midplane, the inward advective flux of vapor is absent in the upper layers of the active disk runs.

To further illustrate this process, we examine the frequency with which particles cross specific radial slices near the snowline region. We select several radial slices, spaced at 0.05 au intervals, that span both sides of r_pk. Each slice is divided into multiple vertical intervals, and we record the number of particle crossings moving towards the star (N₋) and away from the star (N₊). The net crossing count at each slice, defined as N_net = N₊ − N₋, sums to the total number of particles, N_p (1000 in this case).

In Figs. 8c and d, we focus on the snowline region of the active-nqL and passive-nqL runs. The net crossing fraction (N_net/N_p), shown in color, reveals two distinct mechanisms for water transport in the two scenarios. In the active-nqL case, water primarily exits the snowline through the midplane, via diffusion (indicated by the blue blob at 2.0 au), after crossing a barrier formed by the outflow (shown by the red band along the snowline). Conversely, in the passive-nqL run, water particles leave the snowline in addition through upper-layer accretion, as indicated by blue blobs at higher altitudes. The total number of crossing (N_c), indicated by the size of the circles, peak around r_pk (the star symbols), with the active-nqL run exhibiting higher values than the passive-nqL run. This observation is consistent with Figs. 8a,b, indicating that the highly stochastic motion of water parcels near the snowline slows down accretion and that the “upper-layer” channel is significantly more efficient than the “midplane” channel.

6 Discussion

In this section, we discuss the implications of the simulation results on planetesimal formation (Sect. 6.1) and subsequent pebble accretion (Sect. 6.2) in the snowline region. New observational features of snowline derived from the simulation is discussed in Sect. 6.3. Finally, we highlight the limitations of this study and propose potential directions for future research (Sect. 6.4).

Fig. 8 Upper panels: summary of the Lagrangian trajectory integration in the active-nqL and passive-nqL disks. Yellow-green color shading denotes the total time spent in each grid cell by the water particles, which are released from the outer boundary. The black dashed lines denote the root-mean-square of the positions of particles in the vertical direction. The stars denote rpk of different runs (Table 1). The thick grey lines represents the location where ρp,ice/ρvap = 1 and 10−3, respectively, depicting the regions where particles are mostly in ice or vapor phase, respectively (see Appendix A for detailed explanation). The blue and red arrows represent the velocity of ice and vapor, respectively. The length of the arrows is proportional to vHg /D (D, v are the diffusivity and velocity of gas or pebble), respectively, which represent the ratio of the diffusion and advection timescales, $t_{\text{dif}}/t_{\text{adv}}\sim \left(H_{\text{g}}^2/D\right)/\left(H_{\text{g}}/v\right)$. We limit the length of the arrows so that so that vHg/D ≤ 3 for visualization purposes. Also, the depth of color of blue arrows represents tdif,peb/tadv,peb. Lower panels: crossing frequency of Lagrangian particles at various altitudes along the radial direction. Circles are plotted at the center of each vertical interval, with their surface area proportional to the total number of crossings (Nc). Color indicates the ratio of net crossings to the total particle number (Nnet/Np). Here, negative (blue) values denote net accretion, whereas positive (red) values denote net outflow. The domain is zoomed-in to focus on the snowline region. The checkerboard pattern in the net crossing fraction around the snowline seen in the active-nqL run is a manifestation of vigorous water cycling.

6.1 Pile-up and planetesimal formation

Several studies have pointed out that a local solids-to-gas ratio of ∼1 is necessary to trigger SI. The SI would amplify the solids- to-gas ratio by orders of magnitude and result in planetesimal formation through gravitational collapse (Johansen & Youdin 2007; Bai & Stone 2010a; Carrera et al. 2015; Li & Youdin 2021). The water snowline has long been suggested as a promising site to trigger SI due to its ability to enhance the midplane solids- to-gas ratios (Ros & Johansen 2013; Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017; Hyodo et al. 2019). Our study specifically focuses on the “vapor retro-diffusion” scenario, which suggests the formation of icy planetesimals outside snowline, as the “traffic jam”-induced pile-up is minor when pebble disaggregation is omitted (see Fig. 4a and introduction).

Figure 9 summarizes the outcome of the 2D simulations in terms of ξ_pk. Two 1D runs are also included. In 1D/SO17 the gas velocity is fixed to the unperturbed viscous solution and in 1D-hydro the gas motion is solved hydrodynamically (see Sect. 3.2 and Appendix B). For comparison, the estimated solid- to-gas ratio in an unperturbed disk², where no ice sublimation takes place, is also highlighted. Already, in 1D/SO17, the midplane solid-to-gas ratio is elevated at the snowline compared to the unperturbed disk due to the diffusion of vapor. From Fig. 9 it can be seen that the resolved midplane outflow in the 1D-hydro run further elevates ξ_pk by a factor of 1.7, highlighting the importance of the flow pattern.

Comparing the 1D-hydro run with the 2D runs, the latter produce overall stronger pile-ups, as indicated by the larger M_ice and ξ_pk. This highlights the importance of accounting for 2D flow patterns, where preferential accretion of vapor in the upper layers and settling of ice-rich pebbles establishes the water cycle across the snowline region, especially in active disks (Sect. 5.2). Among the 2D runs, the water ice mass stored in the pile-up region exhibits a clear difference (≈1.5 times) between the group of active disks and the group of isothermal and passive disks. This is due to a distinct snowline morphology and resulting distinct strength of the water cycle (see Sect. 5.2). Though the inclusion of latent heat reduces ξ_pk and shifts r_snow,mid inwards, the water mass budget only changes slightly, implying that M_ice is mainly controlled by the snowline morphology.

None of the simulation runs have reached the values of ξ_pk as high as unity to trigger SI. However, inclusion of the back- reaction of solids on the gas and accounting for disaggregation of ice-rich pebble after sublimation could boost the solid-to-gas ratio by factors ≈4 (Schoonenberg & Ormel 2017) or even trigger runaway pile-up (Hyodo et al. 2019). In addition, the headwind velocity is found to be reduced by a factor of ≈0.6 (Sect. 3.4), which would bring down the metallicity threshold to trigger SI perhaps by a similar value (Bai & Stone 2010b; Abod et al. 2019; Baronett et al. 2024). However, this dependence is possibly non-linear and depends on pebble size (Bai & Stone 2010b). Further study is needed to quantify the effect of headwind on SI to provide a complete picture of planetesimal formation at the snowline.

One key result of the hydrodynamical simulation in this work is the discovery of a water cycle, where water parcels are continuously exchanged between the solid and gas phases, enhancing the pile-up of solids (ice). The water-cycle is more pronounced in disks with a large vertical temperature gradient, i.e., active disks with efficient heating in the midplane and high opacity to preserve the heat. Recently, young Class-0/I disks, have been found to have large dust scale height (Villenave et al. 2023; Lin et al. 2023; Guerra-Alvarado et al. 2024b) and systematically higher accretion rate (~10⁻⁷ M_⊙ yr⁻¹). They are therefore likely to host vigorous water cycles. Hyodo et al. (2021) proposes the preferential formation of giant planet’s icy core during the early evolution stage when the disk snowline moves outward into the outer active zone (i.e., away from the inner MRI “dead zone,” e.g., Gammie 1996). They argue that the large diffusivity in the active zone facilitates outward diffusion of vapor and pile-up of ice exterior to the snowline. With the water cycle preferentially operated in young, active disks, our 2D results further support early formation of planetesimals.

Fig. 9 Pile-up properties of 1D and 2D runs as listed in Table 1. For each run, two markers are plotted to indicate the locations of the midplane snowline rsnow,mid (left) and the peak solid-to-gas ratio rpk (right), with their y-axis indicating the value of the peak midplane solid-to-gas ratio ξpk . The vertical dashed line denotes the location of rsnow,mid before pebble injection (2.1 au). The horizontal dotted line, labeled as “unperturbed”, denotes the expected solid-to-gas ratio at the snowline without ice sublimation. The length of the shaded region, centered on rpk , denotes the FWHM of the pile-up. The color represents the total amount of water mass stored, as ice, in the snowline region.

6.2 Onset of pebble accretion

After the formation of planetesimals, it has been suggested that pebble accretion accelerates planetary growth (Ormel & Klahr 2010; Lambrechts & Johansen 2012). Recently, the viability of pebble accretion in planetesimal rings has been explored in scenarios involving either “clumpy rings,” where pebble drift slows down by the back-reaction on the gas (Jiang & Ormel 2023), rings supported by pressure bumps (Morbidelli 2020; Guilera et al. 2020; Lee et al. 2022), or rings in smooth disks (Liu et al. 2019; Jang et al. 2022). To activate pebble accretion, planetesimals need to exceed a threshold mass (Ormel & Klahr 2010; Visser & Ormel 2016; Liu & Ji 2020). Specifically, Liu & Ormel (2018) found that pebble accretion is fully activated at the mass M_* = η³τ_sm_*, which translates to a radius of $$\begin{array}{l}{R}_{*}={\left(\frac{{\tau }_s{\eta }^3{M}_{\star }}{4/3\pi {\rho }_{•,\text{pltsml}}}\right)}^{1/3}\\ \approx 308\text{\hspace{0.17em}km}\left(\frac{\eta }{0.0013}\right){\left(\frac{{\rho }_{•,\text{pltsml}}}{1.0\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-3}}\right)}^{-1/3}{\left(\frac{{\tau }_{\text{s}}}{0.03}\right)}^{1/3}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{1/3}.\end{array}$$(37)

where ρ_•,pltsml is the internal density of planetesimal and η = 0.0013 corresponds to its minimum value at the snowline region in active disks (see Fig. 6). If this mass is not achieved by the planetsimal formation process, the planetesimals first need to grow via collisions in their birth belt (Liu et al. 2019). However, N-body simulations find that the eccentricity and inclination of planetesimals are quickly excited by self-scattering, suppressing the accretion efficiency and making the collisional growth phase the main bottleneck to spawn planets (Liu et al. 2019; Jang et al. 2022; Kaufmann et al. 2025).

Several studies have investigated the size distribution of planetesimals generated from SI (Simon et al. 2016; Schäfer et al. 2017; Abod et al. 2019; Li et al. 2019). Inspired by numerical simulations, Liu et al. (2020) derives a general expression of the characteristic mass of planetesimals, which varies with the stellar mass and aspect ratio. Taking the solar mass and h ≈ 0.04 adopted from our simulation at r_pk, the characteristic radius of planetesimals translates to ≈331 km ≳ R_*,. Therefore, planetesimals formed at the snowline will immediately start efficient pebble accretion. In contrast, a value of η ≈ 0.0023 (Fig. 6) exterior to the snowline region, results in a much larger R_* ≈ 546 km that falls short of the pebble accretion threshold. The reduction of η due to vapor injection (Sect. 3.4) therefore promotes planet formation at the snowline region.

Assuming the 3D limit, at r_pk these planetesimals would grow at a rate ${\dot{m}}_{\text{p}}=12{\tau }_{\text{s}}{m}_{\text{p}}{\rho }_{\text{peb}}{\mathrm{\Omega }}_{\text{K}}{r}^3{f}_{\text{set}}^2/M_{\star }$ (Ormel & Liu 2018), corresponding to a timescale, $$\begin{array}{l}{t}_{\text{PA},3\text{D}}=\frac{m_{\text{p}}}{{\dot{m}}_{\text{p}}}=3.0\times {10}^3{\left(\frac{r}{2.2\text{au}}\right)}^{-3/2}{\left(\frac{{\rho }_{\text{peb}}}{7\times {10}^{-11}\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-3}}\right)}^{-1}\\ \times {\left(\frac{{\tau }_{\text{s}}}{0.03}\right)}^{-1}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{1/2}{\left(\frac{f_{\text{set\hspace{0.17em}}}}{0.62}\right)}^{-2}\text{yr}.\end{array}$$(38)

Here, f_set represents the settling fraction, with f_set ≈ 0.62 for a 331 km-sized planetesimal in our case. The pebble density, ρ_peb = ρ_p,ice + ρ_p,sil, has its value measured from r_pk at the midplane of the active-nqL run. This expression may break down at large m_p if accretion enters the 2D limit or when ṁ_p exceeds the pebble supply from the outer disk (Ṁ_peb). Nevertheless, Equation (38) shows that pebble accretion is rapid, leading to a mass doubling time of only several thousands of years. As a result, we expect the fast growth of a proto-Jupiter’s icy core (≈15 M_⊕) within 0.05 Myr at the snowline region.

6.3 Observational features of the water snowline

Despite the important role of the water snowline in planet formation, pinpointing its exact location within disks remains elusive (e.g., Zhang 2024). Around quiescent stars, it is extremely resolution-demanding to pinpoint the inner 3 au region (e.g., Taurus, 0.02″ for 3 au at 150 pc; Guerra-Alvarado et al. 2024a) where the snowline is thought to reside. Conversely, the out- bursting FU Ori-type (FUor) systems, whose accretion rate can be elevated to ~10⁻⁴ M_⊙ yr⁻¹ (Fischer et al. 2023), will have extended their snowlines to 10–100 au, offering us a unique opportunity to hunt snowlines. Intriguingly, V883 Ori, an FUor, shows an intensity depression (dark annulus) at ≈40 au in ALMA band 6 continuum (Cieza et al. 2016). Initially, this dip was interpreted as a change in dust density and opacity due to accumulation of small grains, which drift slower interior of the snowline. This interpretation requires that the disk is optically thin at millimeter wavelength. On the other hand, Kóspál et al. (2021) argued that FUors are massive, compact, and optically thick. In addition, the depression in line emission from rare isotopes (e.g., HDO, C¹⁷O, Tobin et al. 2023) seen for the inner 40 au of V883 Ori, would imply optically-thick attenuation. Further insight comes from recent multi-wavelength analysis, which reveal that the spectral index remains continuous across the intensity depression, indicating that grain properties are not significantly altered across the dip (Houge et al. 2024). These findings suggest that the snowline region is optically thick and the intensity depression cannot be attributed to variations in dust properties.

These arguments leads to an appealing alternative explanation in line with the findings of this work: the observed intensity depression reflects a dip in temperature due to latent heat cooling. In Fig. 7b, it is observed that the temperature is reduced by ~20% (40 K on top of 200 K) at the sublimation front, compared to model without latent heat effects, and that this temperature plateau affects a region of radial extent of about 50%r_snow,mid. This reduction is sufficient to account for the intensity depression observed in V883 Ori (≈15% in band 7, adopted from Houge et al. 2024).

While it is appealing to connect the intensity dip to latent heat, we caution that outbursting systems bear more nuances than quiescent systems. The accretion luminosity and, hence, the disk temperature, evolves on a relatively short timescale (over about a century in the case of FUors; Fischer et al. 2023), whereas the transport timescale is much longer at the outer region of the disks (Schoonenberg et al. 2017; Colmenares et al. 2024). Therefore, it is likely that the outbursting systems are still evolving (non-steady) at the time of observation. The precise location of the moving snow surface further depends on the sublimation and condensation timescales (e.g., Schoonenberg et al. 2017). From an observational perspective, there is also no consensus on the location of snowline in V883 Ori. Rather than being at 40 AU, molecular line emission suggests that it lies at ≈80 au (Tobin et al. 2023) or even father out, depending on whether the emission layer probes the water sublimation region near midplane (van’t Hoff et al. 2018; Leemker et al. 2021). Recently, Houge et al. (2025) demonstrated that a large amount of water could be hidden by opaque dust from being probed in the infrared band. To probe the temperature plateau, the sublimation region must also extend vertically beyond the emission layer of ALMA continuum. In future work, we would like to include these timedependent aspects and observational uncertainties to understand the observational features of snowline in outbursting systems.

Additionally, the outburst itself could shape dust evolution and sublimation/condensation processes. After the outburst, the size distribution can deviate from collisional equilibrium for a long time (Houge & Krijt 2023). Depending on the cooling efficiency, condensation may occur preferentially on small silicate grains by heterogeneous nucleation rather than on ice-coated pebbles (Ros & Johansen 2024). Future studies should assess the roles of the various physical processes during outburst to obtain a robust interpretation of current observations.

6.4 Future directions

In this work, we have investigated the effect of snowline morphology and latent heat on snowline properties. Previous studies show that the mass accretion rate (Ṁ_acc), ice-to-gas flux ratio (f_i/g) , disk viscosity (diffusivity), pebble size (τ_s) would all influence the pile-up strength, particularly the peak dust-to-gas ratio (Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Hyodo et al. 2019). How these parameters interact with key snowline processes identified in this work, such as the water cycle and latent heat cooling is non-trivial. As the key parameter to drive flow pattern and active heating, we present results for different α in Appendix C, demonstrating the robust existence of the water cycle. While an extensive parameter exploration has been omitted in this work, we discuss their potentially significant influence to guide future study. For example, smaller-sized pebbles are found to enhance pile-ups due to their slower radial drift (e.g., Schoonenberg & Ormel 2017). With the larger scale height and surface area of smaller pebbles, the recapture of vapor could be more efficient. However, the settling of pebbles will be less efficient, leading to uncertainty on the influence of the strength of water cycle when a realistic size distribution of pebbles is considered. Smaller pebbles are also more tightly coupled to the gas, resulting in a significantly reduced or even reversed drift velocity due to the outflow identified near the snowline. In this scenario, silicate grains can be a significant component of the pile-up by the “outflow-induced” traffic jam (see Appendix C), in contrast to the ice-dominated pile-up found in this work (Sect. 3.2). This outflow-induced traffic jam will be further exacerbated if pebble disintegration is considered (e.g., Hyodo et al. 2019). Additionally, coagulation is expected to occur efficiently in regions of solid pile-up (or dust rings), driven by high particle concentration and reduced drift speeds (e.g, Jiang et al. 2024). Understanding how the processes of coagulation and fragmentation shape the dust size distribution will be critical to determining water trapping efficiency near the snowline. Another key parameter is the ice-to-gas flux ratio f_i/g. Typically lower f_i/g leads to smaller ξ_pk of the solid pile-up (e.g., Schoonenberg & Ormel 2017). However, in the active disk, lowering f_i/g also mitigates the latent cooling effect that extends the width and reduces ξ_pk of the pile-up (Sect. 4.1). Consequently, ξ_pk is not necessarily smaller with lower f_i/g . A thorough parameter exploration is necessary to arrive at a quantitative understanding and provide prescriptions for planet formation models relying on a snowline.

The simulation setup in this work can be improved regarding its thermodynamic treatment. First, we assumed that the infrared opacity, which determines the cooling efficiency, is proportional to the non-condensible gas density (H and He) κ = κ₀ρ_H,He (Sect. 2.6). In reality, icy mantle of dust grains could dominate the infrared opacity, given its higher abundance compared to silicates. Accounting for opacity from ice in the context of solar system’s water snowline evolution, Oka et al. (2011) find a slight outward shift of the snowline position due to an enhanced “blanketing effect.” Second, we assumed the standard α-viscosity disk where the heating is deposited in the midplane. However, in MHD simulations focusing on the inner disk region, Mori et al. (2019) finds that the main heating source in the water snowline region is the Joule heating (due to the Ohmic resistivity of disk) that occurs at several gas scale heights above the midplane. In contrast, the pebble scale height in dead zone could be ≪H_g (Zhu et al. 2015; Xu et al. 2017; Yang et al. 2018). We therefore expect that the snowline morphology in MHD disks will resemble the isothermal and passive disks presented in this work, which should be confirmed in future simulations. Since the influence of latent heat on temperature is weaker when heat is not deposited in the midplane (as discussed in Sect. 4.2), we expect that the broadening of solid pile-up and the reduction of ξ_pk by latent heat are also weakened in MHD disks. Moreover, the flow structure in MHD disk is highly complex and different from viscous disk (e.g., Bai 2017), which could drive different levels of pile-up, even becoming asymmetric about the midplane (Hu & Bai 2021).

While steady-state solutions always emerge in our simulations, it was recently proposed that the CO snowline can become thermodynamically unstable when ice sublimation occurs (Owen 2020). This instability arises when stellar heating is optically thick, but midplane cooling is optically thin, causing sublimation to reduce the cooling rate while the heating remains largely unaffected, leading to further temperature increases and more sublimation. The cooling conditions resemble the iso run in this work (Fig. 3, the midplane is optically thin for cooling while κ_vi = 10 cm⁻²g⁻¹ ). However, since we did not account for opacity change due to ice sublimation, the instability cannot be triggered. In addition, a latent heat exchange may operate against the instability given its cooling effect during ice sublimation. Although the latent heat effect was argued to be negligible in Owen (2020), it has been observed to at least locally alter the temperature structure in the iso run (Fig. 3). Therefore, future 2D models incorporating opacity changes and latent heat effects should be conducted to confirm the feasibility of this snowline instability.

7 Conclusions

We conducted 2D multifluid hydrodynamic simulations towards the water snowline. The simulations include pebble dynamics, vapor transport, and phase change processes, while solving the temperature structures self-consistently with a two-stream radiative transfer method. We investigate how different snowline morphology and latent heat effects influence the solid pile-up outcome and explore the possible observational implications. Our main findings are:

1. Viscous spreading of vapor released by inward-drifting ice-rich pebbles at the snowline drives a radial outflow, ≈10 times stronger than the background gas motion (Fig. 5). This previously unrecognized advective outward flux of vapor is of the same strength as the outward diffusive flux, augmenting the pile-up of ice outside the snowline;

2. We identified a water cycle across the snowline region in active disks, where the temperature in the midplane is greater than the overlying layers (Fig. 3). Pebbles sublimate when passing through the snowline region near the midplane. The released vapor diffuses to cooler upper layers, where it recondenses on pebbles and contributes to their growth. These pebbles then settle and return the water back to the midplane (Sect. 3.3). This water cycle reduces the capability of water to escape from the snowline region via the upper layers, which, combined with the midplane outflow, significantly enhances the mass budget of ice in the snowline region (Sect. 5.2). Therefore, active disks with a temperature inversion are more conducive to planetesimal formation than passive or isothermal disks;

3. The total ice mass stored in the snowline region is mainly determined by the strength of the water cycle, which, in turn, depends on the morphology of the snowline (Fig. 9). Radial temperature gradient across the snowline determines the extent and peak of ice pile-up, with a larger width and smaller peak dust-to-gas ratio under flatter temperature profiles;

4. The release of vapor at the snowline reduces the azimuthal velocity of the gas in the disk (η parameter; Fig. 6), which lowers both the metallicity threshold for streaming instability and the threshold radius for planetesimals to accrete pebbles by a factor of ≈0.6–0.7. Planetesimals could immediately start pebble accretion after they are born in the snowline region. Due to the elevated densities, pebble accretion could occur rapidly with optimal growth timescale around 0.05 Myr for giant planets’ icy core (Sect. 6.2);

5. The latent heat exchange tends to flatten the temperature gradient across the snowline, extending the width of the snowline region, while reducing the peak solid-to-gas ratio (Fig. 7). These latent heat effects are most pronounced in active disks with strong blanketing effects;

6. Latent heat cooling causes an extended plateau in the thermal structure of disks, up to 50% of the snowline location (r_snow,mid) in radial extent and ~40 K in temperature (Fig. 7). Such temperature dips would have observational imprints in the dust continuum emission at radio wavelengths (e.g., an intensity dip). In particular, outbursting systems, where the water snowline is pushed to ~10–100 au, are promising targets to observationally test the predictions reported in this work with facilities such as ALMA.

In general, our effort to resolve the snowline in 2D hydrodynamic simulations further supports the snowline as a promising site for the formation of first-generation planetesimals. In the future, incorporating additional physical elements and processes into simulations, along with observational constraints, will enable a more quantitative understanding of the snowline and its significance in planet formation.

Acknowledgements

The authors thank the anonymous referee for providing a thoughtful and helpful report. YW acknowledges helpful discussions with Enrique Macias, Haochang Jiang, Wenrui Xu, Ruobing Dong, Rolf Kuiper, Mario Flock and Anders Johansen. This work is supported by the National Natural Science Foundation of China under grant Nos. 12250610189, 12233004 and 12473065. SM acknowledges support by Japan Society for the Promotion of Science KAKENHI grant no. 22K14081. The authors acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University for providing computational and data storage resources that have contributed to the research results reported within this paper.

Appendix A Lagrangian particle tracking

Our approach to implement the Lagrangian particle tracking follows Ciesla (2010). After a timestep ∆t the positions of water particles are updated through $${(}_r={(}_r+v_{\text{eff}}\mathrm{\Delta }t+R{\left[\frac{2}{\zeta }D\mathrm{\Delta }t\right]}^{1/2},$$(A.1)

with ζ = 1/3, D the diffusivity of water in either ice or vapor state (we drop the subscript “p” or “g” hereafter for simplicity) and R a random number between [–1, 1] following uniform distribution. The effective velocity v_eff = (v_eff,r, v_eff,ɀ) has only its ɀ-component given in Ciesla (2010). Here we re-derive it to obtain also the r-component.

The dynamical evolution of water parcels follows the diffusion-advection equation, $$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v-D{\rho }_{\text{g}}\nabla c\right)=0,$$(A.2)

where c = ρ/ρ_g is the concentration of either ice or vapor. In cylindrical coordinate, the equation reads, $$\frac{\partial \rho }{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\rho \right)+\frac{\partial }{\partial ɀ}\left(v_{ɀ}\rho \right)=\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_{\text{g}}Dr\frac{\partial c}{\partial r}\right)+\frac{\partial }{\partial ɀ}\left({\rho }_{\text{g}}D\frac{\partial c}{\partial ɀ}\right).$$(A.3)

To reproduce the dynamical evolution of water parcels subject to diffusive-advective motions in the way of Monte Carlo tracking of particles, the Monte Carlo methods should satisfy the Fokker- Planck equation: $$\frac{\partial f}{\partial t}+{\displaystyle \sum _{i=1}^N\frac{\partial }{\partial x_i}}\left(v_{i}f\right)={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N\frac{{\partial }^2}{\partial x_i\partial x_j}}}\left(D_{ij}f\right),$$(A.4)

where i, j denote different dimensions and N = 2 in our case. The Fokker-Planck equation governs the evolution of the probability distribution of the variable ρ, which is denoted as f in above equation. In the simulation we assume D_ij = D(r), then the equation becomes: $$\begin{array}{l}\frac{\partial f}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_{r}f\right)+\frac{\partial }{\partial ɀ}\left(v_{ɀ}f\right)=\frac{1}{r}\frac{\partial }{\partial r}\left[r\frac{\partial (Df)}{\partial r}\right]+\frac{{\partial }^2}{\partial {ɀ}^2}(Df)\\ \equiv \frac{1}{r}\frac{\partial }{\partial r}\left[r\frac{\partial D}{\partial r}f+rD\frac{\partial f}{\partial r}\right]+\frac{{\partial }^2}{\partial {ɀ}^2}(Df).\end{array}$$(A.5)

To rewrite the diffusion-advection equation into the form of Fokker-Planck equation, we expand Equation (A.3) as $$\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\rho \right)+\frac{\partial }{\partial ɀ}\left(v_{ɀ}\rho \right)=\frac{1}{r}\frac{\partial }{\partial r}\left(rD\frac{\partial \rho }{\partial r}\right)\\ --\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{D}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial r}\rho \right)+\frac{{\partial }^2}{\partial {ɀ}^2}(D\rho )-\frac{\partial }{\partial ɀ}\left(\frac{D}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial ɀ}\rho \right).\end{array}$$(A.6)

Comparing Eqs. (A.5) and (A.6), we can get $$\begin{array}{l}{v}_{\text{eff},r}=v_r+\frac{D}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial r}+\frac{\partial D}{\partial r},\\ v_{\text{eff},ɀ}=v_{ɀ}+\frac{D}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial ɀ}.\end{array}$$(A.7)

The second terms in both v_eff,r and v_eff,ɀ comes from density stratification, that particles are more likely to appear near denser region. The last term in v_eff,r comes from the radial gradient of diffusivity. In the Monte Carlo simulation, v_r and v_ɀ are taken from the steady-state velocity field of either ice or vapor in the hydrodynamic simulations.

Next, the phase change probability is determined according to the sublimation rate (R_e) and condensation rate (R_c). Following Ros & Johansen (2013), we can define R_e and R_c as $$\frac{dm}{dt}=4\pi R^2{v}_{\text{th\hspace{0.17em}}}{\rho }_{\text{vap\hspace{0.17em}}}\left(1-\frac{P_{\text{eq}}}{P_{\text{vap\hspace{0.17em}}}}\right)\equiv R_{\text{c}}{\rho }_{\text{vap\hspace{0.17em}}}-R_{\text{e}}{\rho }_{\text{p},\text{ice}},$$(A.8)

where m is the mass of a pebble. Therefore, $$\frac{R_{\text{c}}}{R_{\text{e}}}=\frac{P_{\text{vap\hspace{0.17em}}}}{P_{\text{eq}}}\frac{{\rho }_{\text{p,ice\hspace{0.17em}}}}{{\rho }_{\text{vap\hspace{0.17em}}}}.$$(A.9)

Then, the probability of having phase changes taking place for Lagrangian particles follows as $$P=\{\begin{array}{cc}{R}_{\text{e}}/\left(R_{\text{e}}+R_{\text{c}}\right),& (\text{\hspace{0.17em}sublimation\hspace{0.17em}}),\\ R_{\text{c}}/\left(R_{\text{e}}+R_{\text{c}}\right),& \text{\hspace{0.17em}(condensation)}.\end{array}$$(A.10)

In Fig. 8, we plot ρ_p,ice/ρ_vap to denote the respective probability to be either in ice or vapor phase, since R_c/R_e closely follows ρ_p,ice/ρ_vap at the region where phase changes happen frequently.

Appendix B Update of 1D snowline model in gas dynamics

To verify our hydrodynamic model, we simulated the 1D snowline following the setup in Schoonenberg & Ormel (2017). The disk temperature profile is fixed as $$T(r)=150{\left(\frac{r}{3\text{\hspace{0.17em}au}}\right)}^{-1/2}\text{\hspace{0.17em}K}.$$(B.1)

The gas is supplied at the outer boundary at a constant rate Ṁ_acc = 10⁻⁸M_⊙ Yr⁻¹. Other parameters (ƒ_i/g, τ_s,0, α) are also kept the same as 2D runs (Table 1). In this way, the disk snowline is expected to be located at ~2.1 au. In previous studies (Schoonenberg & Ormel 2017; Hyodo et al. 2021), the gas and dust motion are solved with a set of transport equations, where the gas velocity and the surface density of non-condensible gas (Σ_xy) are fixed as the viscous solution of an unperturbed disk composed of pure H-He, $$v_{\text{g}}=\frac{3\nu }{2r};{\mathrm{\Sigma }}_{\text{xy}}=\frac{{\dot{M}}_{\text{acc}}}{3\pi \nu }.$$(B.2)

This treatment ignores two important aspects which would reshape the solid pile-up. Firstly, with ice sublimation, the surface density of the gas is altered with injected vapor and the gas velocity should deviate accordingly following (e.g., Lynden-Bell & Pringle 1974), $$v_{\text{g}}=-\frac{3}{r\mathrm{\Sigma }{v}_{\text{k}}}\frac{\partial }{\partial r}\left(r\mathrm{\Sigma }{v}_{\text{k}}\nu \right).$$(B.3)

Secondly, as vapor diffusion occurs, the H-He gas also diffuses in a similar manner. Consequently, Σ_xy no longer sticks to the pure advection solution. Instead, it must be determined by the combined effects of advection and diffusion, as described by the advection-diffusion process (see Equation (8)).

By incorporating the momentum equation and solving for the diffusion of all species within the hydrodynamic model, the changes in gas velocity and the effects of H-He diffusion are inherently accounted for. To demonstrate their effect on the pile-up, a “simple model” (1D/SO17) is constructed, where both the gas velocity and Σ_xy are fixed following Equation (B.2). In Fig. B.1, we compare the pile-up strength of the “simple” model and the contrasting “full” model (1D-hydro). As expected, the surface density of H-He gas (Σ_xy) deviates from the unperturbed viscous solution, manifesting a density bump outside the snowline (Fig. B.1a). Importantly, the 1D-hydro run presents around a factor of two’s increase in the peak solid-to-gas ratio (Fig. B.1b). These can be understood as the result of midplane viscous outflow of gas near the snowline, which significantly slows down the drift of pebbles (Fig. B.1c) and leads to a stronger pile-up. The 1D results therefore confirm the importance of gas dynamics in shaping the solid pile-up outside the snowline, which is already discussed in Sect. 5.1.

Fig. B.1 Comparisons of 1D/SO17 and ID-hydro. (a): the surface density profiles. The H-He surface density is timed by the ice-to-gas ratio ƒi/g for a clear comparison. (b): the solid-to-gas ratio at the disk midplane and the vertically-integrated surface density ratio. The grey dashed line denotes snowline location (2.1 au). (c): the gas and pebble radial velocity at the disk midplane. It is worth noticing that the gas velocity in 1D/SO17 haven’t returned to the unperturbed viscous solution even at the outer boundary, giving rise to the difference of Σxy all over the place outside the snowline.

Appendix C Effect of disk viscosity

To illustrate the effect of viscosity, we vary the disk viscosity parameter α to 0.03 and 0.001 with respect to the fiducial value 0.003 (we do not choose α = 3 × 10⁻⁴ since it takes too much computational effort to reach a steady state). Here we only study the active disk case, where the water cycle is vigorous. Following the procedure described in Sect. 2.6, the initial snowline is anchored at 2.1 au; therefore, κ_R can be adjusted accordingly. All other parameters (e.g., Ṁ_acc and f_i/g) are kept the same as in the active run (see Sect. 2 and Table 1). The steady state radial profiles are presented in Fig. C.2 and the density contours in Fig. C.1.

First, the water cycle persists with various disk viscosity, as can be seen from the settling of pebbles in Fig. C.1. Second, as α increases, the water cycle extends over a larger region in the r − ɀ plane. This extension is also clearly seen in surface density profiles of ice, where the pile-ups become wider with increasing α (Fig. C.2). This occurs because (1) the increased scale height of pebbles allows vapor to recondense at higher altitudes, and (2) the stronger gas background outflow in the midplane — driven by higher disk viscosity (e.g., Takeuchi & Lin 2002) — impedes pebble drift. Comparing the movement of pebbles in Fig. C.1a,b and Fig. 2b, we observe that despite a fixed Stokes number, pebbles drift more efficiently at lower α, whereas at α = 0.03, they are strongly influenced by gas motion and can even move outward near the midplane. As a result, a higher α leads to much higher level of solid pile-up in terms of both the amount of ice accumulated and the solid-to-gas ratio (Fig. C.2).

In the case of α = 0.03, the strong outflow not only promotes the pile-up of ice, but also leads to the pile-up of silicate. As shown in Fig. C.2 (left panel), the surface density of silicate at the snowline becomes comparable to that of ice. This arises due to a traffic jam effect induced by the gas outflow, which slows silicate drift near the snowline. Consequently, predictions based purely on the vapor retro-diffusion scenario may significantly underestimate the strength of solid pile-up. For instance, in steady state, our simulation yields ξ_pk ≈ 1.5 for α = 0.03, whereas Schoonenberg & Ormel (2017) reports ξ_pk ≈ 0.8 (see their Fig. 8). Future works should take into account the effect of gas hydrodynamics, especially at high viscosity.

Fig. C.1 Steady state density structure of active disk runs of different α. The notation is the same as Fig. 2 except for different α, ρ0 (the midplane gas density at r0 of the unperturbed disk) is different (ρ0 ≈ 3.4 × 10−12 cm g−3 for α = 0.03 and ≈1.0 × 10−10 cm g−3 for α = 0.001).

Fig. C.2 Steady state radial profiles of active disk runs of different α. The notation is the same as Fig. 4. Given the fixed Ṁacc , the gas and surface density is much larger when α = 0.001 since they are accreted slower.

References

All Tables

All Figures

Fig. 1 Initial midplane snowline location resulting from the two-stream radiation transport model (Equation (29)) as function of Rosseland mean opacity κR and stellar luminosity L⋆ . Here the optical opacity is fixed to κvi = 10 cm2g−1. The permitted parameter space is divided into three regimes (white dashed lines, for illustration purpose only): “passive disk” when L⋆ is high, “active disk” when κR is high and “vertically-isothermal” disk when κR is small. The dotted lines are temperature contours of the midplane at 2.1 au, ranging from 100 K to 200 K in 10 K increments. The black solid line represents the parameter combinations that lead to rsnow,mid = 2.1 au and the black dots denote the parameter combinations used in this work, as summarized in Table 1. Hydrodynamic and thermodynamic effects will move the snowline location in the simulation.

In the text

Fig. 2 Steady state density structure of all simulated disk (as listed in Table 1). The dark green line denotes the boundary where ρp,ice/ρg = 10−2. The region enclosed by the green line filled with blue shading represents the ice density normalized to ρ0 ≈ 3.4 × 10−11 cm g−3 (the midplane gas density at r0 of the unperturbed disk), while the region interior to the snowline indicates the vapor density with red shading. The ice (blue) and vapor (pink) mass fluxes are plotted with streamlines, whose thickness denotes the magnitude of the flux normalized to ρ0cs,0. The mass flows in the vertical direction, which are especially vigorous in the active-nqL and active runs, contribute towards an increased pile-up of ice in the snowline region.

In the text

Fig. 3 Steady state temperature structure of all the simulation runs (Table 1). The color shows the temperature structure with contours highlighting the range between 150 and 200 K in increments of 5 K. The dashed lines show the Rosseland mean optical depth τR from 0.1 to 10.0, displaying the optically-thick or optically thin nature of the cooling. In panel (e) and (f), the dashed lines are absent since the disk is very optically thin.

In the text

Fig. 4 Radial profiles of the passive run. (a): vertically-integrated surface density of total gas, ice, silicate and vapor. The grey line shows the expected vapor surface density from the unperturbed viscous solution (Equation (2)), where fi/g is the ice-to-gas flux ratio. The blue shade denotes the region that contains Mice. (b): the solid/vapor-to-gas ratio in the midplane. The positions of peak solid-to-gas ratio (rpk) and the midplane snowline (rsnow,mid) are highlighted; (c): vertically- integrated radial mass fluxes of ice, vapor, water, H-He and silicate, where ℱwater = ℱvap + ℱice. The dashed lines denote the imposed gas (Ṁacc) and ice ( fi/g Ṁacc) fluxes, respectively.

In the text

Fig. 5 Upper panels: vertically-integrated radial mass flux of vapor. For three runs (passive,activeand active-nqL), the diffusive, advective and total vapor flux are plotted. The dotted lines denote zero-flux levels. The blue arrows indicate the position of the midplane snowline (rsnow,mid). The advective flux contributes significantly to the total vapor flux in theactivedisk runs (active and active-nqL), while it is negligible in the passive run. Lower panels: midplane gas radial velocity and the contributions from different terms following Equation (36). The thick grey line represents the measured value from the simulation. The T2 term signifies the dominant contribution of the radial viscous stress to the radial velocity, while the contributions from radial advection (T1 ) and ɀ-direction (Tɀ) are negligible.

In the text

Fig. 6 Azimuthal velocity deviation at the midplane , expressed in terms of the radial pressure gradient parameter η (Equation (35)). The black dashed and solid lines show η before vapor injection and after the steady state is reached, respectively. The red lines show the vapor fraction corresponding to the y-axis on the right. The blue arrows indicate the position of the midplane snowline (rsnow,mid). A reduction of η at the snowline region strengthens solid pileups and is conducive to planetesimal formation by streaming instability and planet formation by pebble accretion.

In the text

Fig. 7 Comparisons of radial profiles of the active disks. For all three panels, the solid line represents the active disk while the dashed line represents active-nqL disk. (a) The surface density profiles. As in Fig. 4, the expected vapor surface density fi/g Ṁacc /(3πν) is plotted in grey. (b) The temperature profile of slices in θ-direction at the midplane and higher altitudes. The slices above the midplane pass through a certain height ɀ at the location of snowline (rsnow) of the active run. (c) The solid-to-gas ratio (black) and vapor-to-gas ratio (red) at the midplane in steady state.

In the text

Fig. 8 Upper panels: summary of the Lagrangian trajectory integration in the active-nqL and passive-nqL disks. Yellow-green color shading denotes the total time spent in each grid cell by the water particles, which are released from the outer boundary. The black dashed lines denote the root-mean-square of the positions of particles in the vertical direction. The stars denote rpk of different runs (Table 1). The thick grey lines represents the location where ρp,ice/ρvap = 1 and 10−3, respectively, depicting the regions where particles are mostly in ice or vapor phase, respectively (see Appendix A for detailed explanation). The blue and red arrows represent the velocity of ice and vapor, respectively. The length of the arrows is proportional to vHg /D (D, v are the diffusivity and velocity of gas or pebble), respectively, which represent the ratio of the diffusion and advection timescales, $t_{\text{dif}}/t_{\text{adv}}\sim \left(H_{\text{g}}^2/D\right)/\left(H_{\text{g}}/v\right)$. We limit the length of the arrows so that so that vHg/D ≤ 3 for visualization purposes. Also, the depth of color of blue arrows represents tdif,peb/tadv,peb. Lower panels: crossing frequency of Lagrangian particles at various altitudes along the radial direction. Circles are plotted at the center of each vertical interval, with their surface area proportional to the total number of crossings (Nc). Color indicates the ratio of net crossings to the total particle number (Nnet/Np). Here, negative (blue) values denote net accretion, whereas positive (red) values denote net outflow. The domain is zoomed-in to focus on the snowline region. The checkerboard pattern in the net crossing fraction around the snowline seen in the active-nqL run is a manifestation of vigorous water cycling.

In the text

Fig. 9 Pile-up properties of 1D and 2D runs as listed in Table 1. For each run, two markers are plotted to indicate the locations of the midplane snowline rsnow,mid (left) and the peak solid-to-gas ratio rpk (right), with their y-axis indicating the value of the peak midplane solid-to-gas ratio ξpk . The vertical dashed line denotes the location of rsnow,mid before pebble injection (2.1 au). The horizontal dotted line, labeled as “unperturbed”, denotes the expected solid-to-gas ratio at the snowline without ice sublimation. The length of the shaded region, centered on rpk , denotes the FWHM of the pile-up. The color represents the total amount of water mass stored, as ice, in the snowline region.

In the text

Fig. B.1 Comparisons of 1D/SO17 and ID-hydro. (a): the surface density profiles. The H-He surface density is timed by the ice-to-gas ratio ƒi/g for a clear comparison. (b): the solid-to-gas ratio at the disk midplane and the vertically-integrated surface density ratio. The grey dashed line denotes snowline location (2.1 au). (c): the gas and pebble radial velocity at the disk midplane. It is worth noticing that the gas velocity in 1D/SO17 haven’t returned to the unperturbed viscous solution even at the outer boundary, giving rise to the difference of Σxy all over the place outside the snowline.

In the text

	Fig. C.1 Steady state density structure of active disk runs of different α. The notation is the same as Fig. 2 except for different α, ρ0 (the midplane gas density at r0 of the unperturbed disk) is different (ρ0 ≈ 3.4 × 10−12 cm g−3 for α = 0.03 and ≈1.0 × 10−10 cm g−3 for α = 0.001).
In the text

	Fig. C.2 Steady state radial profiles of active disk runs of different α. The notation is the same as Fig. 4. Given the fixed Ṁacc , the gas and surface density is much larger when α = 0.001 since they are accreted slower.
In the text