The midinfrared water vapor emission spectrum provides a novel way to characterize the delivery of icy pebbles toward the innermost lt5 au regions of planetforming disks. Recently JWST MIRIMRS showed that compact disks exhibit an excess of lowenergy water vapor emission relative to extended multigapped disks suggesting that icy pebble drift is more efficient in the former. We carry out detailed emissionline modeling to retrieve the excitation conditions of rotational water vapor emission in a sample of four compact and three extended disks within the JWST Disk Infrared Spectral Chemistry Survey. We present twotemperature HSUB2SUBO slab model retrievals and for the first time constrain the spatial distribution of water vapor by fitting parametric radial temperature and column density profiles. Such models statistically outperform the twotemperature slab fits. We find a correlation between the observable hot water vapor mass and stellar mass accretion rate as well as an anticorrelation between cold water vapor mass and submillimeter dust disk radius confirming previously reported water line flux trends. We find that the midIR spectrum traces HSUB2SUBO with temperatures down to 180300 K but the coldest 150170 K gas remains undetected. Furthermore the HSUB2SUBO temperature profiles are generally steeper and cooler than the expected superheated dust temperature in passive irradiated disks. The column density profiles are used to estimate icy pebble mass fluxes which suggest that compact and extended disks may produce markedly distinct innerdisk exoplanet populations if local feeding mechanisms dominate their assembly.
['Protoplanetary disks', 'Planet formation', 'Water vapor', 'James Webb Space Telescope', 'Astrochemistry', 'Infrared spectroscopy', '1300', '1241', '1791', '2291', '75', '2285', 'Astrophysics - Earth and Planetary Astrophysics']
Retrieval of Thermally Resolved Water Vapor Distributions in Disks Observed with JWSTMIRI
{'Retrieval of Thermally-Resolved Water Vapor Distributions in Disks Observed with JWST-MIRI': "<!-- image -->  \n9  \n1 Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA 02138, USA 2 Department of Physics, Texas State University, 749 N Comanche Street, San Marcos, TX 78666, USA 3 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr. Pasadena, CA, 91109, USA 4 Vassar College, 124 Raymond Avenue, Poughkeepsie, NY 12604, USA 5 NSFs NOIRLab, 950 N. Cherry Avenue, Tucson, AZ 85719, USA 6 Division of Geological & Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA 7 School of Physics and Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK 8 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey, RH5 6NT, UK 10 Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA ∗  \n11 Dipartimento di Fisica, Universit'a degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy", 'ABSTRACT': "The mid-infrared water vapor emission spectrum provides a novel way to characterize the delivery of icy pebbles towards the innermost ( < 5 au) regions of planet-forming disks. Recently, JWST MIRI-MRS showed that compact disks exhibit an excess of low-energy water vapor emission relative to extended multi-gapped disks, suggesting that icy pebble drift is more efficient in the former. We carry out detailed emission line modeling to retrieve the excitation conditions of rotational water vapor emission in a sample of four compact and three extended disks within the JDISC Survey. We present two-temperature H 2 O slab model retrievals and, for the first time, constrain the spatial distribution of water vapor by fitting parametric radial temperature and column density profiles. Such models statistically outperform the two-temperature slab fits. We find a correlation between the observable hot water vapor mass and stellar mass accretion rate, as well as an anti-correlation between cold water vapor mass and sub-mm dust disk radius, confirming previously reported water line flux trends. We find that the mid-IR spectrum traces H 2 O with temperatures down to 180-300 K, but the coldest 150-170 K gas remains undetected. Furthermore the H 2 O temperature profiles are generally steeper and cooler than the expected 'super-heated' dust temperature in passive irradiated disks. The column density profiles are used to estimate icy pebble mass fluxes, which suggest that compact and extended disks may produce markedly distinct inner-disk exoplanet populations if local feeding mechanisms dominate their assembly.", '1. INTRODUCTION': 'Ever since the first detection of water vapor in a protoplanetary disk by Carr et al. (2004) (see also Carr & Najita 2008; Salyk et al. 2008), there has been considerable interest in using its infrared emission spectrum as a tracer of the thermochemical and dynamical properties  \nCorresponding author: Carlos E. Romero-Mirza  \ncarlos.romero mirza@cfa.harvard.edu  \n- ∗ NASA Hubble Fellowship Program Sagan Fellow  \nof the innermost ( < 5 au) disk regions. One particularly interesting trend that emerged from the Spitzer InfraRed Spectrograph (Spitzer-IRS) observations is an anti-correlation between sub-mm dust disk masses and mid-IR H 2 O line luminosities relative to those of organics (Carr & Najita 2011; Najita et al. 2013). This was interpreted as evidence that more massive disks are more efficient at forming icy planetesimals, which would end up locking up most of the O-rich ice. Using resolved continuum imaging with ALMA, a follow-up study by Banzatti et al. (2020) in turn indicated that the strongest  \nanti-correlation occurs between mid-IR H 2 O line fluxes and the sub-mm pebble disk radius-and this correlation persists even after accounting for the accretion luminosity of the host. Hydro-dynamical processes have provided a physical explanation for such a trend: if, unlike multi-gapped extended disks, compact disks are created by efficient pebble drift (e.g. Pinilla et al. 2012; Appelgren et al. 2020; Toci et al. 2021; van der Marel & Mulders 2021), brighter H 2 O line emission could indeed be observed as due to the increased rate of ice sublimation near the snowline (Cyr et al. 1998; Ciesla & Cuzzi 2006; Kalyaan et al. 2021, 2023).  \nA more direct observational test to determine whether compact disks exhibit efficient pebble drift can be obtained by measuring the excess emission of moderateenergy (with upper-level energies E up ≲ 3000 K) rotational water lines. This is because such H 2 O transitions are expected to be in near local thermo-dynamic equilibrium (LTE), and can be dominated by emission from gas temperatures between 200-400 K (Salyk et al. 2011; Du et al. 2017; Banzatti et al. 2020, 2023a) consistent with the expected thermal conditions in elevated layers near the snowline, while warm ( > 400 K) gas-phase water formation would dominate the emission from higher energy levels (Glassgold et al. 2009; Meijerink et al. 2009; Blevins et al. 2016; Bosman et al. 2022). Measuring the relative excitation of lower- vs higher-energy water lines was not possible earlier with Spitzer-IRS spectra, since most mid-IR H 2 Olines were blended together at its low spectral resolving power ( R ∼ 600). Instead, with JWST MIRI-MRS (henceforth MIRI-MRS) (Wells et al. 2015; Argyriou et al. 2023) it is now possible to carry out spectroscopic studies of water vapor with unprecedented detail thanks to its sensitivity, wavelength coverage, and resolving power (2000 -4000) (e.g. Grant et al. 2023; Gasman et al. 2023; Pontoppidan et al. 2024; Xie et al. 2023). Indeed, MIRI-MRS can now be used to search for direct evidence as to whether compact disks exhibit brighter low-energy H 2 O emission lines compared to extended disks (Banzatti et al. 2023b).  \nThe aim of this work is to go beyond comparing the H 2 O spectral line patterns between compact and extended disks Banzatti et al. (2023b), and carry out detailed emission line modeling of the entire rotational spectra to retrieve the excitation conditions of water vapor in a sample of compact and extended disks with a number of dust gaps. We present a novel technique to fit H 2 O emission spectra using not only single- (e.g. Grant et al. 2023) or two-temperature (Pontoppidan et al. 2024; Romero-Mirza et al. 2024) slabs, but also using physically-motivated parametric models to characterize the radial temperature and column density struc-  \nture of water vapor in each disk. It is known that single slab models cannot reproduce the entire rotational H 2 O emission spectrum of disks due to excitation temperature ( T ex ) gradients (Blevins et al. 2016; Liu et al. 2019; Banzatti et al. 2023a; Gasman et al. 2023), but so far radial gradients have rarely been included in midIR molecular excitation retrievals (Zhang et al. 2013; Kaeufer et al. 2024). By modeling temperature and density gradients, we can start to use H 2 O emission as a proxy to study the inner disk physical structure-from otherwise spatially- and kinematically-unresolved observations. Extracting radial distributions, furthermore, makes it possible to provide estimates for the inner-disk icy pebble mass fluxes, which in turn can give clues as to which kind of planets can form via the pebble accretion mechanism across localized radii in the disk (Ida et al. 2016; Lambrechts et al. 2019).  \nThis paper is organized as follows: in Section 2, we present the MIRI-MRS data used in this paper. Section 3 presents the slab modeling technique and Bayesian retrieval details, as well as their results. We discuss our findings in Section 4 and, finally, present our conclusions in Section 5.', '2. OBSERVATIONS': 'We analyze new and published MIRI-MRS observations of seven protoplanetary disks around T Tauri stars. Four systems have compact pebble disks (with resolved sub-mm dust disk radii R dust < 30 au) while three are significantly more extended ( R dust > 100 au). Our sample includes the four H 2 O-rich disks first presented by Banzatti et al. (2023b): GK Tau, CI Tau, IQ Tau, and HP Tau, which provided the first confirmed cases of cold water vapor excess in compact disks. In addition, we analyze the compact disk of FZ Tau, known to harbor an abundant reservoir of hot and cold water vapor (Pontoppidan et al. 2024), and the large disk of AS 209, which exhibits bright dust continuum emission and a mid-IR spectrum dominated by hot H 2 O Romero-Mirza et al. (2024). Finally, we present new MIRI-MRS observations of the compact disk around the classical T Tauri star (cTTs) GQ Lup.  \nThe disk of GQ Lup was observed as part of JWST Cycle 1 program GO-1640 (P.I. Andrea Banzatti) on August 13, 2023, with deep integrations of 1700 -1800 s per MRS sub-band, with the specific goal of characterizing the degree of inner-disk water vapor enrichment. For observational details regarding the rest of our sample we direct the reader to Banzatti et al. (2023b), Pontoppidan et al. (2024), and Romero-Mirza et al. (2024). The MIRI-MRS data of all disks in our sample were processed using the JWST Disk Infrared Spectral Chem-  \nistry Survey (JDISCS) reduction pipeline version 7.1, which includes an empirical de-fringing process using asteroid calibrators. This is fully described by Pontoppidan et al. (2024). Table 1 summarizes the properties of the disks in our sample.  \nFigure 1 presents the continuum-subtracted spectra for all the disks included in this study, and the wavelength region considered in our analysis. The dust continuum is subtracted by fitting a Gaussian Process model to the data, as described in Romero-Mirza et al. (2024). We also experiment with using a Savitzky-Golay filter (e.g. Pontoppidan et al. 2024) to subtract the continuum emission, and find negligible differences. A plot showing full spectra and empirical continuum models for each disk is presented in Appendix A.  \nIn this paper, we focus on the analysis of water vapor rotational emission from each of these disks. The identification and retrieval of column densities and excitation conditions of organics and other species in these and other sources within the JDISCS survey will be presented in an upcoming overview paper (Arulanantham et al., in prep.).', '3. MODELING AND RESULTS': 'We aim to retrieve the physical properties of water vapor emitting from the innermost regions of the seven disks in our sample. The mid-IR emission spectrum of H 2 O is generated using the LTE slab modeling Python package iris (Romero-Mirza 2023) with spectroscopic data obtained from the HITRAN database (Gordon et al. 2022). The slab models consist of opacity-weighted intensities that account for line blending and line-center saturation effects, which is critical to obtain accurate column density measurements. The opacity-weighted H 2 O line intensities are modeled assuming a line Gaussian width given by the sum in quadrature of a 1 km/s turbulent component and a thermal broadening component √ k B T ex /µ , where µ is the molecular weight of H 2 O. Typical thermal widths are of 0.4-0.8 km/s for the best-fit temperatures reported below. We only consider in our modeling the 12 . 0 -27 . 0 µ m wavelength range, since emission at shorter wavelengths is dominated by high-energy ro-vibrational H 2 O lines likely emitting from hot and diffuse disk regions where water vapor is not thermalized by collisions with H 2 (Carr et al. 2004; Meijerink et al. 2009; van Dishoeck et al. 2014; Banzatti et al. 2023a). iris has been benchmarked against the publicly available package spectools ir (Salyk 2022).  \nThe slab models are calculated assuming a different resolving power within the wavelength coverage of each MIRI-MRS sub-band. First, we generate the line emission models at a sufficiently high resolution ( R ≈ 10 5 ) to sample the line profiles. Then, the models are convolved with a Gaussian kernel of width ∆ λ = λ m /R , where λ m is the median wavelength of each sub-band, and binned onto a wavelength grid with uniform spacing ∆ λ / 2. To calculate the log-likelihood of the models, the MIRI-MRS data and their uncertainty in each subband is binned onto the same uniform wavelength grid while ensuring flux conservation using the Python package SpectRes (Carnall 2017). Appendix B describes how the flux uncertainty in each sub-band is estimated. Table 2 lists the resolving power and wavelength sampling assumed for each instrument sub-band. We emphasize that the full 12 . 0 -27 . 0 µ m region is modeled simultaneously.  \nThe modeling and retrieval details are presented below. In brief, we first model the H 2 O emission in each disk as the sum of two slabs to represent two temperature components: a first approximation to a radial temperature gradient as done e.g. in Pontoppidan et al. (2024). Next, we fit power law excitation temperature and column density profiles to each disk. Finally, based on the results of the two previous retrievals, we fit a more flexible model with a tapered power law column density profile. All models are fit using the Bayesian nested sampling Python package dynesty (Speagle 2020). We apply multiple ellipsoid decomposition and random walk sampling, and use the change in remaining evidence as a measure of convergence.', '3.1. Two-Temperature Water Vapor Models': "Observations with Spitzer-IRS , Herschel-PACS , and ground observatories first provided indications that the IR spectrum of water vapor from disks cannot be typically reproduced by a slab model with a single excitation temperature, but rather multiple temperature components are necessary (for detailed overviews see van Dishoeck et al. 2014; Banzatti et al. 2023a). This has recently been confirmed with slab model fits to MIRIMRS disk spectra (Banzatti et al. 2023b; Gasman et al. 2023; Xie et al. 2023; Pontoppidan et al. 2024), finding that models with two temperature components can reproduce the rotational H 2 O line emission reasonably well. We first model the H 2 O spectra for the disks in our sample following this simpler approximation.  \nSpecifically, the spectra are modeled as the sum of two plane-parallel slabs, each with an emitting area ( A ), water column density ( N ), and excitation temperature  \nFlux + shift (Jy)  \nFigure 1. Continuum-subtracted MIRI-MRS spectra of each of the disks in the sample, vertically shifted for clarity. The blue shaded regions indicate the wavelength regions excluded from our analysis due to ro-vibrational emission ( < 12 µ m) or larger noise ( > 27 µ m).  \n<!-- image -->  \n( T ex ) (Romero-Mirza et al. 2024). We use the following Uniform (in log 10 space) prior distributions for A and N :  \nlog 10 N = U (15 , 20) cm -2 log 10 A = U ( -3 , 3) au 2 , (1)  \nand the following Normal priors for the excitation temperatures:  \nlog T = N (2 . 9 , 0 . 1) K (2)  \nlog 10 T cold = N (2 . 6 , 0 . 2) K 10 hot  \nTable 1. Properties of the disk sample  \nNote -Distances obtained from GAIA DR3 Gaia Collaboration et al. (2023). Disk inclinations ( i ) obtained from MacGregor et al. (2017); Long et al. (2019); Pontoppidan et al. (2024). Sub-mm dust disk radii ( R dust ) obtained from Huang et al. (2018); Long et al. (2019, 2020); Wu et al. (2017), and correspond to the radius enclosing 95% of the emission. Stellar luminosities ( L ⋆ ) and masses ( M ⋆ ) obtained from Andrews et al. (2018); Herczeg & Hillenbrand (2014); Alcal'a et al. (2017); McClure (2019). Mass accretion rates obtained from Salyk et al. (2013); Herczeg & Hillenbrand (2014); Alcal'a et al. (2017); McClure (2019). Innermost dust gap locations from Huang et al. (2018); Long et al. (2018, 2020). a Putative cavity based on the near infrared index (see discussion in Banzatti et al. 2023b). b The inner gap found in GQ Lup is very shallow and may not produce an effective trap against pebble drift (Kalyaan et al. 2023). c An additional inner gap at 5 au has been proposed by Jennings et al. (2022) using super-resolution techniques.  \nTable 2. Slab model details  \nNote -The resolving power and wavelength sampling used to model the water vapor emission in each wavelength region of interest. The values of R are estimated from the empirical fits reported by Pontoppidan et al. (2024).  \nfor the cold and hot water vapor components, respectively. These correspond to Normal priors centered around 400 and 800 K, with standard deviation of ∼ 200 K. Note that the use of Normal priors does not affect the best-fit parameters in any meaningful way, and they are simply used to consistently guide one slab towards the low- and another towards the high-temperature solution. The priors are informed by previous best-fit water vapor parameters reported by Banzatti et al. (2023b), Pontoppidan et al. (2024), and Romero-Mirza et al. (2024) for six of the disks analyzed in this work.  \nFigure 2 presents the best-fit emitting areas, column densities, and excitation temperatures for each disk. Here and in the rest of this paper, the 'best-fit' values are chosen to be the median of the posterior distribution obtained for each parameter. Table 3 lists the exact values and corresponding uncertainties for each parameter. We find that while the emitting areas of the cold and hot water vapor reservoirs are different from each other, they are each similar across the sample. The cold H 2 O components tend to have emitting areas of 1 -10 au 2 (or equivalent radii R eq = 0 . 6 -1 . 8 au assuming circular areas), whereas the hot components have areas of order 0 . 1 au 2 ( R eq = 0 . 2 au). The hot H 2 O column densities are also fairly similar for compact and extended disks, with typical values of 10 18 -10 19 cm -2 . On average, however, we find that compact disks exhibit higher cold H 2 O column densities than extended disks, by a factor of a few. With respect to the excitation temperatures, we retrieve comparable T ex values across the sample of 800 -1000 K for hot H 2 O and ∼ 400 K for cold H 2 O, with some variation: HP Tau and CI Tau have warmer cold H 2 O components ( ∼ 500 K), while AS 209 exhibits a very cold H 2 O component ( ∼ 250 K).  \nNote that the poor quality of the AS 209 spectrum affects the robustness of this particular fit. In agreement with Romero-Mirza et al. (2024), we find that there are multiple high-energy H 2 O lines with sufficient SNR to accurately constrain the hot water vapor component. Yet, while detected, the lower-energy H 2 O lines are af-  \nFigure 2. Best-fit two-temperature H 2 O model parameters for each disk in our sample. Error bars correspond to the 0.16 and 0.84 percentiles. Observable masses are calculated as the product of column density, emitting area, and molecular weight. The disks are listed in ascending R dust order from left to right, and the dashed lines separate compact and extended disks . Note that sometimes the error bars are smaller than the marker.  \n<!-- image -->  \nfected by noise artifacts beyond 20 µ m. This results in a significantly higher relative uncertainty for the retrieved emitting area and column density of cold water vapor in AS 209 compared to the rest of the sample (see Figure 2). The excitation temperature, however, appears to be well constrained by the line ratios. As further discussed in Section 4.4, while the peculiar properties of cold H 2 O in AS 209 may arise due to the lower SNR, they could  \nalso indicate the presence of inner disk substructure, for which there is kinematic evidence (Banzatti et al. 2023a).  \nGiven the variation in excitation temperatures among hot and cold components in the sample, we caution against over-interpreting apparent correlations (or lack thereof) between H 2 O column densities and/or areas and disk sizes at this stage. That is, some of the water vapor classified as 'cold' in some disks would contribute to the 'hot' component in others. While one could in principle alleviate this by fitting the data with fixedtemperature components (e.g. one with T ex = 400 and one with T ex = 800 K), we find that doing so produces significantly larger residuals. This, combined with the fact that we do not observe strong correlations between T ex , A , and N , seems to indicate that the source-tosource variation in retrieved temperatures occurs due to distinct radial temperature structures between disks.  \nFigure 3 shows the best-fit two-temperature H 2 O model for the compact disk of GQ Lup. The contributions of the hot and cold components are shown individually, and the spectrum has been broken up into multiple panels for clarity. We emphasize that including a cold H 2 O component is critical to reproduce the low-energy water lines beyond 20 µ m. Note that most of the bright lines not reproduced by the model (e.g. near 16.8, 18.8, 20.1, 21.5, 23.2, and 25.3 µ m) are OH emission lines. By contrast, Figure 4 presents the bestfit two-temperature H 2 O model to the extended disk of CI Tau. Note how in this case, as reported in Banzatti et al. (2023b), the spectrum is dominated by hot H 2 O even at long wavelengths, with little excess from colder H 2 O emission. The two-temperature models fit to the rest of the sample are presented in Appendix C.", '3.2. Temperature and Column Density Gradients': 'It is likely that the two-temperature H 2 O fits are tracing the average physical conditions representing radial excitation temperature and abundance gradients in each disk, which are expected to contribute to the rotational water spectrum given the large spread of upperlevel energy levels spanned by the lines in MIRI spectra (Banzatti et al. 2023b). Thus, we explore whether there is enough information in the MIRI-MRS spectra to retrieve the properties of these radial gradients, and whether straightforward parametric models can provide improved model fits to the data relative to using twotemperature models. The details are presented below.  \n3.2.1. Power law Profiles  \nTable 3. Two-component H 2 O fit results  \nNote -The reported values for each parameter consist of the median and 0.160.84 posterior distribution percentiles.  \nTo fit a radial distribution of water vapor, the first necessary assumption we make is that the excitation temperature profiles can be described by a monotonically decreasing function of radius. We choose to approximate the underlying H 2 O excitation temperature distribution using a power law (e.g. Kaeufer et al. 2024) of the form  \nT ex = T 0 ( r 0 . 5 au ) -α . (3)  \nNote that if the H 2 O gas is thermalized with the dust, one would expect a ∼ 1 / 2 for a passive irradiated disk (Chiang & Goldreich 1997; Dullemond et al. 2001). For simplicity, we further assume that the distribution of water vapor is smooth and also described by a power law:  \nN = N 0 ( r 0 . 5 au ) -β , (4)  \nwhich could either increase or decrease with radius depending on dynamical, chemical, and/or optical depth effects.  \nTo implement these parametric models, we generate the H 2 O spectra as the sum of 50 slabs representing concentric rings spanning from 0.1 to 10 au in each disk. The minimum radius of 0.1 au is an estimate based on the typical dust sublimation radius for the disks in our sample. We decide to construct the radial grid using uniformly log 10 -spaced points, in order to more finely sample the innermost 1 au region. Each ring has an emitting area given by:  \nA = π ( r 2 out -r 2 in ) cos i, (5)  \nwhere i is the disk inclination, and r out and r in are the outer and inner ring radii, respectively. The exci-  \ntation temperature and column density of each ring is evaluated at the mid-point.  \nFor these retrievals, we use the following Uniform prior distributions:  \nlog 10 T 0 = U (2 . 0 , 3 . 3) K log 10 N 0 = U (12 , 22) cm -2 α = U (0 . 0 , 5 . 0) β = U ( -5 . 0 , 5 . 0) , (6)  \nto allow the models a reasonably large amount of flexibility. Note that while these models consist of 50 slabs, they require fewer free parameters than the twotemperature fits (only four instead of six).  \nThe models achieve convergence for all disks in the sample, with statistically well-constrained parameters. Figure 5 shows the best-fit power law column density and excitation temperature profiles for each disk. The plotted profiles are truncated at the radius where T ex drops below 150 K, which we identify as the predicted surface snowline location. We find that extended disks exhibit steep column density profiles, with β = 1 . 8 -2 . 1. On the other hand, compact disks exhibit a wide variety of column density profiles: the best-fit profiles for GK Tau and GQ Lup are similar, with β = 0 . 9 -1 . 0, while that of HP Tau is much steeper ( β = 1 . 8), akin to those of extended disks. Interestingly, FZ Tau exhibits a peculiar, radially-increasing H 2 O column density profile, with β = -1 . 1, that peaks at ∼ 5 × 10 19 cm -2 prior to freeze-out. As for the excitation temperature profiles, all disks exhibit profiles with α = 0 . 5 -1 . 0, steeper than expected for an optically-thin, radiationdominated disk layer. We find an average a of 0.73,  \nFigure 3. Best-fit two-temperature H 2 O model to the compact disk of GQ Lup. The data is shown in black, the total model in violet, and the individual contributions of the cold and hot H 2 O components are shown in blue and magenta, respectively, shifted vertically for clarity. Regions where the emission is blended with OH or organics are labeled.  \n<!-- image -->  \nwith no clear temperature profile distinctions between compact and extended disks.', '3.2.2. Exponentially-tapered Column Density': "One potential drawback of fitting power law column density profiles is that in order to reproduce the de-  \npletion of cold water vapor in extended disks-and the sharp depletion of water vapor at the snowline in any disk-the profiles may converge to steeper power law indices ( b ), which in turn increase the column density of hot water vapor. To explore whether the column densi-  \nFigure 4. Same as Figure 3, for the extended disk of CI Tau. Note that unlike GQ Lup, the spectrum of CI Tau is dominated by the hot H 2 O component even at longer wavelengths.  \n<!-- image -->  \nties of high-temperature H 2 O found for extended disks are only artificially higher due to a restrictive parametric model, we next include an exponential taper to the column density profiles:  \nN = N 0 ( r 0 . 5 au ) -β exp [ -( r r taper ) ϕ ] , (7)  \nwhere the taper radius r taper is a free parameter. Implementing the additional parameter ϕ thus allows us to  \nFigure 5. Power law H 2 O column density (top) and excitation temperature (bottom) profiles fit to each disk in the sample. The solid lines correspond to the median profiles generated from posterior samples, while the shaded regions indicate the 0.16 and 0.84 percentiles at each radius. The profiles are truncated at 150 K. The dashed-line portions correspond to the regions where the modeled H 2 O emission lies below the detection threshold (see Section 4.1).  \n<!-- image -->  \nmodel the depletion of cold water vapor without necessarily forcing the column density profile to increase in the innermost disk regions. For these fits, we maintain the Uniform priors described in Section 3.2.1 and use the following Uniform priors for r taper and ϕ :  \nlog 10 r taper = U ( -1 . 0 , 1 . 0) , ϕ = U (0 . 0 , 5 . 0) . (8)  \nThe best-fit profiles are presented in Figure 6. Appendix F shows example corner plots illustrating the posterior distributions obtained for FZ Tau. We find that including an exponential taper helps model the suppression of cold H 2 O in extended disks while preventing the hot H 2 O column densities from reaching unlikely high values ( > 10 19 cm -2 ). The N profiles of three compact disks, GK Tau, HP Tau, and GQ Lup converge  \nFigure 6. Same as Figure 5, for the exponentially-tapered column density profile fits.  \n<!-- image -->  \nto combinations of ϕ and r taper such that they only vary minimally compared to the previous simple power law profile fits. However, the column density profile of FZ Tau changed significantly. The profile now increases as N ∝ r 0 . 5 out to a radius of around 0.7 au, beyond which N drops steeply with ϕ ≈ 2 . 4. Note that in this case adding a taper has also prevented N from reaching extremely high values out to the snowline. Finally, we do not observe significant changes to the excitation temperature profiles, except for that of FZ Tau, which is now less steep with a down from 0.8 to 0.6.  \nIn Figure 7, we present a comparison between the two-temperature, power law, and exponentially-tapered power law best-fit models to FZ Tau. Note that all three prescriptions provide excellent fits to the data, and most differences are only noticed by close inspection of the residuals. For instance some line fluxes are underpredicted by the two-temperature model (see the region near 23.8 µ m) and others are over-predicted (see 18.55, 19.6, 20.35 µ m), while the tapered-power law model can better reproduce these lines. Distinctions like these are observed across the sample, but the emission lines mentioned above are not always affected in the same way.  \nµ  \nFZ Tau  \nWavelength (  \nFigure 8 shows the maximum log-likelihood achieved with each model and the range of log L spanned by the 0 . 025 -0 . 975 (2 σ ) percentiles of posterior samples for  \n<!-- image -->  \nm)  \nFigure 7. Best-fit H 2 O emission models to the spectrum of the compact disk of FZ Tau. The two-temperature model is shown in cyan (bottom), the power law model in blue (middle), and the exponentially-tapered power law model in violet (top of each panel). The data is shown three times in each panel and the models are shifted vertically for clarity.  \nSimilar comparisons are presented for each disk in Appendix D.  \nFigure 8. The maximum log likelihood (circles) attained for each H 2 O emission model fit. The labels 2T, PL, and TPL correspond to the two-temperature, power law, and exponentially-tapered power law models, respectively. The vertical bars show the range of log L values spanned by posterior distribution samples within the 0.025-0.975 percentiles. Note that fitting a temperature gradient always results in a statistically significant improvement, except in the case of AS 209.  \n<!-- image -->  \nevery disk in the sample. The percentiles provide one metric to evaluate which models provide better fits to each dataset. For FZ Tau, we note that a power law N profile provides a worse fit compared to the twotemperature model, but a tapered power law then results in an improved fit. In general we find that fitting T ex and N gradients results in better model fits relative to simply using two-temperature models, with > 95% confidence. This result is in agreement with Kaeufer et al. (2024), who also find statistically significant evidence for an H 2 O column density and temperature gradient in the disk of GW Lup. Yet the specific parametric model used (power law versus tapered power law) only results in statistically significant fit improvements for FZ Tau and CI Tau. Still, the exponentially-tapered profiles are preferred over the power law profiles given that the latter occasionally produce unlikely high column density values 1 . AS 209 is the only case in which both gradient fits perform significantly worse than the two-temperature model. This is further discussed in Section 4. Note that the comparison presented in Figure 8 is rather conservative, as both models with radial gradients not only achieve higher likelihoods, but also require the same or fewer free parameters. The detailed treatment of T ex and N gradients is therefore preferred to accurately model water emission, especially in cases when the H 2 O spectrum must be subtracted to identify and fit blended, less-abundant species.  \nAn additional benefit of fitting N and T ex gradients is that we can now develop a well-defined measure of 'cold'  \nand 'hot' water vapor in each disks. In particular, we define cold water vapor as that with 150 ≤ T ex < 400 K, and hot water vapor that with T ex ≥ 400 K. The motivation behind the 400 K choice is twofold: First, 400 K accounts for radial and vertical diffusion of water vapor from the sublimation front. Second, 400 K is approximately 100 K above the temperature at which the gas-phase formation of water starts to become efficient (Meijerink et al. 2009), thus ensuring that most of the water vapor classified as 'hot' originates from gas-phase chemistry rather than sublimation. The observable cold and hot water vapor masses can then be found by integrating the column density profiles within the corresponding radial regions. While the cold components may be contaminated by gas-phase formation of water, and the hot components by sublimation, we experimented using temperature cut-offs of up to 500 K and find no significant differences in the trends presented below.  \nFigure 9 shows the derived observable cold and hot H 2 O masses for each disk, and as a function of mass accretion rate and outer sub-mm dust disk radius. These masses were calculated using the exponentially-tapered column density profiles. A Spearman correlation test reveals a statistically significant anti-correlation between cold water vapor mass and R dust , and a tentative correlation between hot water vapor mass and the stellar mass accretion rate. The former exhibits a correlation coefficient ρ = -0 . 89 with p -value 0.01, and the latter has ρ = 0 . 75 with p = 0 . 05. In fact, the observed relations can be well-fit by power law expressions of the  \nFigure 9. Left : Observable cold (150 < T ex < 400 K) and hot ( T ex > 400 K) water vapor mass in each disk, calculated by integrating the best-fit exponentially-tapered power law column density profiles over each temperature range. Middle: Observable H 2 O masses as a function of stellar mass accretion rates. The Spearman correlation coefficient and p-value for the hot water mass is indicated on the top. The dashed line indicates a power law fit. Right : Same as middle panel, but for H 2 O masses vs. outer sub-mm dust disk radius. The correlation coefficient corresponds to the cold H 2 O mass.  \n<!-- image -->  \n/circledot  \nTable 4. Radial H 2 O profile fit results  \nNote -The reported values for each parameter consist of the median and 0.16-0.84 posterior distribution percentiles.  \nform  \nlog 10 ( M [cold H 2 O] g ) = p log 10 ( R dust au ) + q, log 10 ( M [hot H 2 O] g ) = r log 10 ( ˙ M M ⊙ /yr ) + s, (9)  \nwith p = -1 . 23 ± 0 . 2, q = 23 . 8 ± 0 . 4, r = 0 . 48 ± 0 . 1, and s = 25 . 6 ± 0 . 9 as shown by the dashed lines in Figure 9. We consider the latter correlation tentative given the larger p-value, and because the high accretion rate of FZ Tau has a strong effect in driving the correlation. No significant correlation is found between cold water vapor mass and mass accretion rate, nor between hot water vapor mass and R dust , even after dividing out the correlation with ˙ M . Finally, to complement Figure 9, Figure 26 in Appendix E shows the average column density ( ˆ N ) and total emitting area for cold and hot water vapor in each disk.", '4.1. Is the Snowline Detected?': 'We have carried out detailed retrievals to estimate the excitation conditions of water vapor in the innermost regions of seven protoplanetary disks, including the effects of radial gradients in excitation temperature and column density. As mentioned in Section 3.2.1, the profiles are truncated at the location where T ex drops below 150 K, approximately the location of the expected surface H 2 O snowline. However, we further examine whether or not H 2 O emission down to 150 K is in fact detected by MIRI in each disk. To do this, we calculate the total wavelength-integrated flux contributed by each of the 50 slab models fit to each disk, for both the power law and tapered power law profiles. We then compare these to the estimated background noise (see Appendix B), and find the radial/temperature bin at which the modeled integrated H 2 O emission drops below a 5 σ detection level.  \nWe denote these as emission cutoff radii ( R out ) and cutoff excitation temperature ( T out ). We interpret these radii as the outermost mid-IR H 2 O emission location in each disk.  \nFigure 10 shows the best-fit exponentially-tapered power law H 2 O model for the disk of GK Tau, with the individual contributions of each slab spectrum shown individually. The models above the detection threshold are shown in green, while those undetected are shown in orange. In this case, R out = 1 . 7 au, and emission down to ∼ 180 K is detected, which lies near the expected range for H 2 O freeze-out. We detect H 2 O emission with T out down to ∼ 200 K in two other disks: GQ Lup and IQ Tau. For the rest, the mid-IR cutoff occurs at temperatures ≳ 250 K. The precise mid-IR cutoff radius and temperature for each disk are reported in Table 5. Note that the profiles shown in Figures 5 and 6 also indicate the location of R out in each disk.  \nOur findings indicate that there are multiple H 2 O emission lines within the MIRI-MRS wavelength coverage that are sensitive to gas with temperatures < 200 K. However, these lines are weak, and are only detected in a few particularly water-rich disks. In the rest of the sample, the abundance of observable water vapor with T ex < 200 K is too low to be detected in the mid-IR. Thus, the reported location of the 150 K surface H 2 O snowline in each disk should be considered a prediction, and its properties must be confirmed in the future by far-IR observations sensitive to the lowest-energy H 2 O lines.  \nWe emphasize that truncating the column density profiles at either 150 K or T out has no significant effect on the cold and hot water vapor mass-radius correlations reported in Section 3.2.2. Hence, whether a) the bestfit N profiles can be extrapolated to T ex < T out or b) the H 2 O column density in this region drops due to e.g. chemical quenching (Blevins et al. 2016) or vertical diffusion (Meijerink et al. 2009) does not affect our conclusions.  \nTo complement Figure 10, Figure 11 shows the integrated H 2 O fluxes as a function of radius for each disk, as well as R out , T out , and the predicted location of the surface snowline. The interpretation of R out is currently open. One possibility is that the mid-IR emission cutoff radius is an indicator of the midplane snowline, which could be the case under the vertical cold-finger scenario (e.g. Meijerink et al. 2009; Krijt et al. 2016) and thus the ratio of R out to the T ex = 150 K radius may provide a measure of vertical transport of gas in a disk. However, the derived R out are inconsistent with the expected location of the midplane snowlines based on viscous heating (Min et al. 2011; Mulders et al. 2015), and are further  \nuncorrelated to ˙ M , contrary to what one would expect (Garaud & Lin 2007).  \nRather, we find a notable distinction between the R out of compact and extended disks. As presented in the bottom panel of Figure 11, the mean H 2 O cutoff radii for extended disks lies at around 0.8 au, whereas that of compact disks is located closer to 2 au. The vertical axis shows the excitation temperature reached at R out , illustrating that, on average, most of the H 2 O flux in compact disks originates from gas with T ex ≳ 220 K, while in extended disks most of the flux comes from water vapor with T ex ∼ 50 K warmer. R out thus appears to be tracing the zone where gas-phase formation of H 2 O becomes efficient at ≳ 300 K in extended disks, whereas in compact disks R out may trace the radius where the column of H 2 O sublimated from icy pebbles becomes detectable in the mid-IR (but not necessarily the theoretical location of the midplane snowline).', '4.2. Comparison with the Dust Temperature': "The derived T ex profiles can be compared to the expected temperature structure in each disk. In particular, we compare the retrieved H 2 O temperature profiles to a simple prescription for the dust temperature above the photosphere-from where the inner-disk water vapor spectrum is expected to arise (e.g. Woitke et al. 2018; Bosman et al. 2022):  \nT dust = 1 ϵ 1 / 4 ( R ∗ 2 r ) 1 / 2 T ∗ , (10)  \nwhere R ∗ and T ∗ are the stellar radii and effective temperature, and ϵ is the dust emissivity in the mid-IR (Chiang & Goldreich 1997; Dullemond et al. 2001). Following Pontoppidan et al. (2024), we assume ϵ = 0 . 17 for the dust size distribution of a standard T Tauri disk, and also explore the effects of setting ϵ = 1 . 0 for dust growth beyond a few µ m and ϵ = 0 . 058 for an ISM dust distribution.  \nFigure 12 shows the modeled dust temperature profiles, compared to the H 2 O excitation temperature profiles fit to each disk. Two results stand out: First, the standard dust temperature appears to be always moderately higher than the H 2 O excitation temperature. Second, the H 2 O T ex profiles are steeper than the dust temperature profiles. The only disk in which the retrieved H 2 O T ex is in very good agreement with the dust temperature is FZ Tau, which is in line with the comparison reported by Pontoppidan et al. (2024).  \nMost thermochemical models predict that the gas and dust temperatures are decoupled in the uppermost disk layers. However, when decoupled, PAH photo-electric heating is predicted to raise the gas temperature by hundreds to thousands of K (e.g. Kamp & Dullemond 2004;  \nFigure 10. The contributions to the H 2 O spectrum fit to the GK Tau data at each radial bin, for the exponentially-tapered column density model. The excitation temperature of each component is indicated on the right side. Models in green have integrated fluxes > 5 σ , whereas the models in orange lie below the detection threshold. In this case, we detect H 2 O emission with T ex down to ∼ 180 K.  \n<!-- image -->  \nWavelength (  \nµ  \nm)  \nJonkheid et al. 2004; Woitke et al. 2009). Evidently, we find no indication of H 2 O gas super-heating. Yet it is unclear what mechanism would lead the gas to be colder than the dust. Assuming that the gas responsible for rotational H 2 O emission is thermalized with H 2 and T ex ≈ T gas , our results may rather simply indicate that gas and dust are coupled and Equation 10 is overpredicting the standard dust temperature. As suggested by Figure 12, the dust temperatures could be lower and  \nin better agreement with the typical H 2 O T ex if there had been significant dust growth in the inner regions of these disks. The steepness of some profiles may further arise due to a radial gradient in dust sizes. Alternatively, the dust temperature could also be lower than expected due to geometric effects such as shadowing from the inner rim (Dullemond et al. 2001; Rafikov & De Colle 2006). Ultimately the MIRI continuum emission may be used for detailed modeling of the dust properties and  \nFigure 11. Top: Normalized wavelength-integrated H 2 O flux profiles as a function of radius, extracted from the bestfit tapered N profiles. Asterisks indicate the radii where T ex reaches 150 K in each disk. Circles show the cutoff radius at which the integrated emission drops below 5 σ . Bottom: The H 2 O cutoff radius for each disk, and the H 2 O excitation temperature derived at that radius. Dashed lines indicate the average values for compact (blue) and extended (red) disks.  \n<!-- image -->  \nout  \ntemperature in each of these disks, but such analysis is beyond the scope of this paper.  \nThe retrieved T ex profiles are generally well constrained, with power law indices a ≳ 1 / 2, steeper than expected for the dust based on radiative equilibrium alone. Interestingly, the average temperature power law index of 0 . 69 found for our sample is in excellent agreement with that predicted by thermochemical models ( ∼ 0 . 7) including the effects of H 2 O self-shielding and chemical heating from photo-dissociation and molecule formation (Bosman et al. 2022). In this context, the H 2 O gas could be cooler than the dust by a few hundred K if H 2 O self-shielding is highly efficient and chemical heating processes are suppressed. Note that the average power law index is also in good agreement with that  \nTable 5. Outer H 2 O cut-off parameters  \nNote -From left to right: The predicted surface snowline location based on the exponentiallytapered models, the outer emission cutt-off radius, and the corresponding H 2 O excitation temperature at that radius.  \nfound for the surface temperature of the GW Lup disk (index of 0.63), based on an H 2 O temperature gradient fit to the MIRI spectrum (Kaeufer et al. 2024).  \nViscous heating may also play some role in setting the H 2 O temperature profiles. While heating from accretion is generally expected to be deposited near the midplane, most of the viscous heat could rather be released near the H 2 O-emitting disk surface via the magnetorotational instability (Hirose & Turner 2011; ' Ad'amkovics et al. 2014). A simple power law may then be insufficient to reproduce the underlying temperature profile, and a more flexible model may be necessary to account for viscous and radiative heating-dominated regions.  \nFinally, it is possible that viscous heating and dust growth/disk shadowing effects collude in such a way that the H 2 O T ex profiles are in fact tracing the dust temperature. That is, the H 2 O emitting region is dominated by a mix of viscous and radiative heating, and the dust temperature modeled by Equation 10 may be an overestimate. This would explain the possible correlation between hot water vapor and mass accretion rates, the steep T ex profiles, and why FZ Tau-which has the highest ˙ M -is the only disk where H 2 O T ex and Equation 10 are in good agreement.", '4.3. The Ring-Like Water Vapor Distribution of FZ Tau': 'The best-fit tapered power law N radial profile reveals an uncommon H 2 O column density distribution in the  \nFigure 12. The excitation temperature profiles (blue) derived for H 2 O in each disk with 1 σ confidence regions, compared to the expected disk surface dust temperature (magenta) for different dust size distributions (see Section 4.2).  \n<!-- image -->  \ndisk of FZ Tau. As opposed to those from the rest of the sample, the N profile in FZ Tau exhibits a ring-like distribution, with N ∝ r 1 / 2 out to 0.7 au and an abrupt exponential drop-off that brings N below 10 16 cm -2 well within the expected freeze-out location (based both on our fits and Herschel observations (Blevins et al. 2016)). It appears unlikely that this phenomenon is somehow the result of inner-disk dust substructure (i.e. a cavity), given the particularly high log 10 ˙ M = -6 . 5 of FZ Tau and its abundant cold water vapor reservoir.  \nConversely, it may be that the high accretion rate of FZ Tau leads to the observed ring-like H 2 O distribution. If the midplane temperature exhibits a large gradient and the radial velocity of icy pebbles is large enough, the bulk of the latter might drift past the location of the water snowline and sublimate in a noticeably warmer disk region (e.g. Piso et al. 2015). This could lead to the the paucity of water vapor with T ex < 300 K, yet still produce a significant enrichment in H 2 O with T ex > 400 K. The small drop in H 2 O column density closer to the star could then be explained by the very high temperatures produced by viscous heating leading to the thermal dissociation of some water vapor (or direct photo-dissociation from the radiation field produced by the accretion itself). Indeed, kinematically-resolved ro-vibrational lines have been observed at 5 µ m with iSHELL from the ground (Banzatti et al. 2023a), and found to have a FWHM of 40 km/s (corresponding to a Keplerian orbit of 0.2 au) and a temperature of 1000 K,  \nmatching well the hotter part of the temperature gradient estimated in this work. Furthermore, Salyk et al. (2019) also find evidence of a ring-like water vapor distribution in FZ Tau based on TEXES observations of 12 µ m H 2 O line emission. In this case the rotational H 2 O lines are consistent with an inner emission cut-off radius of 0.4 au, where we also find evidence of a drop in H 2 O column density.', '4.4. The Temperature Structure of AS 209': 'AS 209 is the disk where we observe the strongest discrepancy between the H 2 O excitation temperature profile and the expected dust temperature from radiative equilibrium. The retrieved T ex profile is the steepest in our sample, and the excitation temperatures are always several hundred K lower than T dust . As discussed above, the steep power law index closer to 1 may suggest that the gas temperature is dominated by viscous heating. However, this is difficult to reconcile with the very low temperatures, even if T dust from Equation 10 was over-estimated.  \nNote that AS 209 is the only disk for which the twotemperature H 2 O model fit statistically outperforms both T ex gradient fits, and the cold H 2 O component has a much lower temperature ( ∼ 250 K) than those found in the rest of the sample. Thus, a more speculative possibility is that the water vapor structure in AS 209 is indeed not well described by a continuous radial distribution: A gas gap somewhere within the innermost 1 au could explain why the disk of AS 209 appears to be de-  \neted in water vapor with T ex within 400 -600 K, but still exhibits a typical hot (800 K) water vapor reservoir and cool 250-300 K component that could be tracing the region beyond the proposed gap. Recent high-resolution iSHELL observations to ro-vibrational CO in AS 209 indeed support the presence of an inner gap from the significant discontinuity in velocity between a broad and narrow component of CO (Figure 18 in Banzatti et al. (2023a)). Using the line velocity, we can estimate that a narrow gap would be somewhere between 0.16 au (from the peak of the broad component) and 0.56 au (from the width of the narrow component). The narrow CO component is also cold (possibly as cold as the cold H 2 O in MIRI), as it significantly drops already at J ∼ 25 while in other disks it is usually strong up to higher J levels. Still, the AS 209 MIRI spectrum has the lowest SNR and line-to-continuum ratio of those included in this study, and it is entirely possible the model mismatch is a sideeffect of the low data quality.', '4.5. Icy Pebble Drift': "We have followed-up on the analysis carried out in Banzatti et al. (2023b), and confirmed that the trend between low-energy H 2 O line fluxes and R dust corresponds to a significant anti-correlation between the latter and the observable water vapor mass with 150 < T ex < 400 K. This new result supports the proposed scenario where compact disks have more efficient icy pebble drift, and thus present an increased flux of H 2 O sublimation near the snowline. By contrast the hot water vapor mass does not correlate with the dust disk size, confirming previous results based on line fluxes alone (Banzatti et al. 2023b).  \nAssuming that the cold water reservoir in each disk is dominated by H 2 O sublimation from ices, we can obtain a simple estimate of the pebble mass flux into the inner disk from the observable water vapor mass as  \n˙ M peb ≈ ( M [Cold H 2 O] t H 2 O )( η f ice f H 2 O ) , (11)  \nwhere t H 2 O is the approximate time that sublimated H 2 O resides in the gas phase, η is a correction factor to account for the unobservable water column under the optically thick layer, f ice is the ice fraction of drifting pebbles, and f H 2 O is the fraction of molecular ice in water. Based on thermochemical modeling, η is expected to be of order 10 3 for water vapor with temperatures of 150-400 K (Meijerink et al. 2009; Bosman et al. 2022). We set η = 10 3 . For f ice and f H 2 O we assume 0.20 and 0.80, which are typically assumed values (e.g. Johansen et al. 2021). Finally, we assume t H 2 O is set by the combination of diffusion ( t diff ) and chemical t chem timescales.  \nTable 6. Pebble flux estimates  \nNote -Pebble fluxes estimated from Equation 11. Disk ages obtained from Neuhauser et al. (2005); Andrews et al. (2018); Long et al. (2019); McClure (2019).  \nFor the former, we can estimate the time it takes for water vapor sublimated from ice in the midplane to reach a typical emitting region at one scale height,  \nt diff = H αc s , (12)  \nwhere c s is the sound speed (e.g. Semenov & Wiebe 2011)  \nFor a typical midplane snowline at 1 au and gas temperature of 200 K, t diff can range between 15-150 yr depending on the assumed inner disk viscosity ( α = 0 . 01 -0 . 001). Once in an elevated disk region, thermochemical models suggest water vapor will only reside in the gas phase for a t chem of 10-50 yr (Calahan, J. priv. comm.) before being photo-dissociated. We hence proceed with an assumed t H 2 O of 100 yr.  \nTable 6 provides the assumed disk ages and estimates of the pebble mass flux and total pebble mass deposited in the inner disk ( M peb ). While extended disks have pebble fluxes of a few tens of earth masses per Myr, those of compact disks are an order of magnitude larger. The compact disk of HP Tau appears to be an exception, with a noticeably lower cold water vapor mass and diminished ˙ M peb of 70 M ⊕ /Myr . As discussed in Banzatti et al. (2023b), the near infrared index of HP Tau shows tentative evidence of a dust cavity within the innermost 2 au. As such, HP Tau may have a different physical structure where molecular gas is affected by different irradiation conditions within a dust cavity, and thus exhibits a reduced cold H 2 O reservoir that mimics those of extended disks.  \nThe disk of GQ Lup presents another interesting case. Despite being the largest of the 'compact' disks, it exhibits a pebble mass flux comparable to that of the smallest disk in the sample, FZ Tau. As shown in the rightmost panel in Figure 9, its cold water vapor mass appears to be higher than expected from the R dust anticorrelation inferred from the rest of the sample. Furthermore, its R out is a factor ∼ 2 greater than those of the other compact disks. These anomalies may be explained by the binary nature of GQ Lup, which is known to have a close brown-dwarf-mass companion at ∼ 110 au (Neuhauser et al. 2005): both dynamical simulations (Zagaria et al. 2021a) and observations (Zagaria et al. 2021b) suggest that radial dust drift is substantially faster in binary systems due to tidal interactions, which would explain the higher abundance of cold water and yet a regular hot H 2 O reservoir. The binary interactions may also lead to higher gas temperatures than expected for an isolated disk (Muller & Kley 2012; Picogna & Marzari 2013), which could lead to the particularly extended outer H 2 O radii found for GQ Lup.  \nThe pebble mass flux estimates hint that compact and extended disks not only have distinct inner disk compositions, but may also result in the formation of different exoplanet populations. For instance, Lambrechts et al. (2019) predict that disks may form either small terrestrial or super-Earth planets depending whether the total deposited pebble mass towards the inner disk is ≲ 100 or ≳ 200 M ⊕ , respectively. In this framework, the compact disks in our sample may be forming systems of water-rich super-Earths, while the large disks may only form small rocky planets (see also discussion in van der Marel & Mulders 2021). Note that there is substantial uncertainty regarding the parameters in Equation 11 and hence on the exact pebble flux values. Thus, the key conclusion is that our analysis provides initial evidence that compact disks experience approximately an order of magnitude higher pebble fluxes compared to extended disks. Very high resolution sub-mm observations (e.g. Guerra-Alvarado et al. 2024), together with source-specific thermochemical modeling will be necessary to obtain more accurate inner-disk pebble flux constraints, and refine the connection between the latter and the observable cold H 2 O reservoir.", '5. SUMMARY': "The main results of this work are summarized below:  \n- 1. We fit two-temperature water vapor slab models to the MIRI-MRS spectra of four compact and three extended disks around T Tauri stars. In general, the emission is very well-characterized by the mix-  \n- ture of an extended (1-10 au 2 ) 400 K component and a compact (0.1 au 2 ) 800 K component.\n- 2. We infer the spatial distribution of water vapor by fitting parametric radial excitation temperature and column density profiles. These models achieve even smaller residuals compared to the two-temperature approximation. Thus, when possible, we recommend this kind of modeling to accurately subtract water vapor emission and search for less abundant species in mid-IR disk spectra.\n- 3. We find that the observable cold water vapor mass is anti-correlated with the sub-mm dust disk size, confirming previously reported line flux correlations. The hot water vapor mass is instead tentatively correlated with stellar mass accretion rate, which may set the mid-IR emitting area.\n- 4. The water vapor excitation temperature profiles appear to be steeper and colder than expected for the 'super-heated' dust temperature above the disk photosphere for a passive irradiated disk. The average power law index of 0.7 is in agreement with thermochemical model predictions.\n- 5. The disk of AS 209 is the only one in our sample that can be better modeled by two independent temperature components rather than a continuous H 2 O distribution. This could be due to inner disk substructure.\n- 6. The disk of FZ Tau exhibits a ring-like column density radial profile, which may be explained by its steep radial temperature gradient and high mass accretion rate.\n- 7. Our models indicate that the bulk of the water vapor emits from a smaller and warmer radial region in extended disks compared to compact disks. This difference in emitting regions is not correlated to the mass accretion rate, and supports the icy pebble drift hypothesis.\n- 8. We estimate icy pebble mass fluxes from the observable cold water masses. Compact disks have typical fluxes an order of magnitude higher than those of extended disks. This may suggest compact disks are relatively more likely to form systems of super-earths.", 'ACKNOWLEDGMENTS': 'The JWST data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The  \nspecific observations analyzed can be accessed via doi: 10.17909/q60w-vr09.  \nK.I. O. acknowledges support from the Simons Foundation (SCOL #321183), an award from the Simons Foundation (#321183FY19), and an NSF AAG Grant (#1907653).  \nA portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technol-  \nogy, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).', 'Facilities: JWST (MIRI-MRS)': 'Software: iris (Romero-Mirza 2023), tinyGP (Foreman-Mackey et al. 2024), dynesty (Speagle 2020), astropy (Astropy Collaboration et al. 2013, 2018, 2022), JAX (Bradbury et al. 2018), matplotlib (Caswell et al. 2024)', 'A. FULL SPECTRA': 'Figure 13 presents the full MIRI-MRS spectra and best-fit continuum models. The continuum emission in each disk is fit following the technique described in Romero-Mirza et al. (2024).', 'B. GAUSSIAN NOISE ESTIMATION': 'We use a Gaussian Process to estimate the amplitude of the background Normally-distributed noise in the continuumsubtracted spectra. In particular, we model the data in each of the six wavelength regions used in our analysis using a squared exponential covariance function k ( λ ) with fixed amplitude 10.0 mJy (e.g. Romero-Mirza et al. 2024) and a correlation length scale given by 2 × λ m /R (see Table 2). We then condition this model on the data, and fit for the variance σ 2 . That is, we assume any spectral features with correlation length smaller than the estimated resolving power in each region arise due to Gaussian noise. Inevitably, this results in some real spectral features being classified as background noise, and thus a rather conservative estimate of the flux uncertainty. Figure 14 illustrates the resulting decomposition in the case of HP Tau. Finally, we add a constant 10.0 mJy uncertainty in quadrature to the best-fit σ , which accounts for the typical uncertainty in the continuum baseline fit. This sum is the uncertainty used to calculate the Gaussian log-likelihood in the Bayesian retrievals described in Section 3. The noise estimates are listed in Table 7.  \nTable 7. Gaussian noise  \nNote -The best-fit Gaussian noise (1 σ ) to each wavelength region included in this study.', 'C. TWO-TEMPERATURE FITS': 'In Figures 15-19 we present the best-fit two-temperature H 2 O models for the four disks not shown in Section 3, showing the individual contributions from hot and cold water vapor.  \nFlux (Jy)  \nFZ Tau  \nFigure 13. The full MIRI-MRS spectrum of each disk (black), and the corresponding continuum emission fit (magenta).  \n<!-- image -->', 'D. MODEL COMPARISON': 'In Figures 20-25 we present the comparison between two-temperature, power law, and tapered power law models for the five disks not shown in Section 3.', 'E. COLUMN DENSITY AND EMITTING AREA CORRELATIONS': 'In Figure 26 we show the average column density and emitting area for hot and cold water vapor plotted against R dust and ˙ M . We note that neither parameter presents a statistically significant correlation with either the mass accretion rate or sub-mm dust disk radius, so the correlations reported for the observable mass appear to arise rather from the interplay of both parameters.  \nHP Tau  \nFigure 14. Gaussian Process decomposition of the spectrum of HP Tau. The black line shows the MIRI-MRS continuumsubtracted spectrum, the magenta line shows the conditioned squared exponential kernel, and the blue shaded region the best-fit ± 1 σ region. The latter includes the uncertainty associated with the continuum baseline fit.  \n<!-- image -->  \n24.5', 'F. EXAMPLE POSTERIOR DISTRIBUTIONS': 'In Figures 27, 28, and 29 we present corner plots with the posterior distributions obtained for each of the three fits (two-temperature, power law, and tapered power law) carried out for the disk of FZ Tau.  \nFor the two-component and power-law posterior distributions, we do not observe strong correlations. We note, however, that often the tapered power-law posterior distributions exhibit correlations between parameters, namely those of the tapered column density profiles. This is because multiple combinations of parameters result in almost identical profiles (as shown in Figure 6, the column density profiles themselves are very well constrained). These  \nFZ Tau  \nFigure 15. Same as Figure 3, for the compact disk of FZ Tau.  \n<!-- image -->  \ncorrelations thus do not affect the results and conclusions presented in this paper, but only the interpretability of the some of the column density profile parameters.', 'REFERENCES': 'Flux + shift (Jy)  \nFigure 16. Same as Figure 3, for the compact disk of GK Tau.  \n<!-- image -->  \nGK Tau  \nAndrews, S. M., Huang, J., P\'erez, L. M., et al. 2018, ApJL, 869, L41, doi: 10.3847/2041-8213/aaf741  \nAppelgren, J., Lambrechts, M., & Johansen, A. 2020, A&A, 638, A156, doi: 10.1051/0004-6361/202037650  \nArgyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111, doi: 10.1051/0004-6361/202346489  \nAstropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33,  \ndoi: 10.1051/0004-6361/201322068  \nAstropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74  \nWavelength (  \nm)  \nFigure 17. Same as Figure 3, for the compact disk of HP Tau.  \n<!-- image -->  \nBanzatti, A., Pascucci, I., Bosman, A. D., et al. 2020, ApJ, 903, 124, doi: 10.3847/1538-4357/abbc1a  \nBanzatti, A., Pontoppidan, K. M., P\'ere Ch\'avez, J., et al.  \n2023a, AJ, 165, 72, doi: 10.3847/1538-3881/aca80b  \nBanzatti, A., Pontoppidan, K. M., Carr, J. S., et al. 2023b, ApJL, 957, L22, doi: 10.3847/2041-8213/acf5ec  \nBlevins, S. M., Pontoppidan, K. M., Banzatti, A., et al. 2016, ApJ, 818, 22, doi: 10.3847/0004-637X/818/1/22 Bosman, A. D., Bergin, E. A., Calahan, J., & Duval, S. E. 2022, ApJL, 930, L26, doi: 10.3847/2041-8213/ac66ce Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs, 0.3.13. http://github.com/google/jax  \nIQ Tau  \nFigure 18. Same as Figure 3, for the extended disk of IQ Tau.  \n<!-- image -->  \nCarnall, A. C. 2017, arXiv e-prints, arXiv:1705.05165,  \ndoi: 10.48550/arXiv.1705.05165  \nCarr, J. S., & Najita, J. R. 2008, Science, 319, 1504,  \ndoi: 10.1126/science.1153807  \n-. 2011, ApJ, 733, 102, doi: 10.1088/0004-637X/733/2/102  \nCarr, J. S., Tokunaga, A. T., & Najita, J. 2004, ApJ, 603, 213, doi: 10.1086/381356  \nCaswell, T. A., Sales de Andrade, E., Lee, A., et al. 2024, matplotlib/matplotlib: REL: v3.7.5, v3.7.5, Zenodo,  \ndoi: 10.5281/zenodo.10669804  \nChiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368,  \ndoi: 10.1086/304869  \nCiesla, F. J., & Cuzzi, J. N. 2006, Icarus, 181, 178,  \ndoi: 10.1016/j.icarus.2005.11.009  \nAS 209  \nWavelength (  \nm)  \nFigure 19. Same as Figure 3, for the extended disk of AS 209.  \n<!-- image -->  \nCyr, K. E., Sears, W. D., & Lunine, J. I. 1998, Icarus, 135, 537, doi: 10.1006/icar.1998.5959  \nDu, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98, doi: 10.3847/1538-4357/aa70ee  \nDullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957, doi: 10.1086/323057  \nForeman-Mackey, D., Yu, W., Yadav, S., et al. 2024, dfm/tinygp: The tiniest of Gaussian Process libraries, v0.3.0, Zenodo, doi: 10.5281/zenodo.10463641  \nGaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940 Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606,  \ndoi: 10.1086/509041  \nGK Tau  \nFigure 20. Same as Figure 7, for GK Tau.  \n<!-- image -->  \nGasman, D., van Dishoeck, E. F., Grant, S. L., et al. 2023, A&A, 679, A117, doi: 10.1051/0004-6361/202347005  \nGlassgold, A. E., Meijerink, R., & Najita, J. R. 2009, ApJ,  \n701, 142, doi: 10.1088/0004-637X/701/1/142  \nGordon, I. E., Rothman, L. S., Hargreaves, R. J., et al.  \n2022, JQSRT, 277, 107949,  \ndoi: 10.1016/j.jqsrt.2021.107949  \nHP Tau  \nFigure 21. Same as Figure 7, for HP Tau.  \n<!-- image -->  \nGrant, S. L., van Dishoeck, E. F., Tabone, B., et al. 2023, ApJL, 947, L6, doi: 10.3847/2041-8213/acc44b  \nGuerra-Alvarado, O. M., Carrasco-Gonz\'alez, C., Mac\'ıas, E., et al. 2024, arXiv e-prints, arXiv:2404.04164, doi: 10.48550/arXiv.2404.04164  \nHerczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97, doi: 10.1088/0004-637X/786/2/97  \nHirose, S., & Turner, N. J. 2011, ApJL, 732, L30,  \ndoi: 10.1088/2041-8205/732/2/L30  \nFlux + shift (Jy)  \nFigure 22. Same as Figure 7, for GQ Lup.  \n<!-- image -->  \nGQ Lup  \nHuang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42, doi: 10.3847/2041-8213/aaf740  \nIda, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72, doi: 10.1051/0004-6361/201628099  \nJennings, J., Tazzari, M., Clarke, C. J., Booth, R. A., & Rosotti, G. P. 2022, MNRAS, 514, 6053,  \ndoi: 10.1093/mnras/stac1770  \nJohansen, A., Ronnet, T., Bizzarro, M., et al. 2021, Science Advances, 7, eabc0444, doi: 10.1126/sciadv.abc0444  \nIQ Tau  \nFigure 23. Same as Figure 7, for IQ Tau.  \n<!-- image -->  \nJonkheid, B., Faas, F. G. A., van Zadelhoff, G. J., & van  \nDishoeck, E. F. 2004, A&A, 428, 511,  \ndoi: 10.1051/0004-6361:20048013  \nKaeufer, T., Min, M., Woitke, P., Kamp, I., & Arabhavi,  \nA. M. 2024, A&A, 687, A209,  \ndoi: 10.1051/0004-6361/202449936  \nWavelength (  \nFigure 24. Same as Figure 7, for AS 209.  \n<!-- image -->  \nµ  \nm)  \nKalyaan, A., Pinilla, P., Krijt, S., Mulders, G. D., &  \nBanzatti, A. 2021, ApJ, 921, 84,  \ndoi: 10.3847/1538-4357/ac1e96  \nKalyaan, A., Pinilla, P., Krijt, S., et al. 2023, ApJ, 954, 66,  \ndoi: 10.3847/1538-4357/ace535  \nKamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991,  \ndoi: 10.1086/424703  \nCI Tau  \nFigure 25. Same as Figure 7, for CI Tau.  \n<!-- image -->  \nKrijt, S., Ciesla, F. J., & Bergin, E. A. 2016, ApJ, 833, 285, doi: 10.3847/1538-4357/833/2/285  \nLambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83, doi: 10.1051/0004-6361/201834229  \nLiu, Y., Pascucci, I., & Henning, T. 2019, A&A, 623, A106, doi: 10.1051/0004-6361/201834418  \nLong, D. E., Zhang, K., Teague, R., & Bergin, E. A. 2020, ApJL, 895, L46, doi: 10.3847/2041-8213/ab94a8  \nFigure 26. Top: the average column density for cold (blue) and hot (orange) water vapor, extracted from the tapered column density profiles, as a function of mass accretion rate (left) and sub-mm dust disk radius (right). Bottom: same as top panel, for the total emitting area of hot and cold water vapor.  \n<!-- image -->  \nLong, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1 Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49, doi: 10.3847/1538-4357/ab2d2d MacGregor, M. A., Wilner, D. J., Czekala, I., et al. 2017, ApJ, 835, 17, doi: 10.3847/1538-4357/835/1/17 McClure, M. K. 2019, A&A, 632, A32, doi: 10.1051/0004-6361/201834361 Meijerink, R., Pontoppidan, K. M., Blake, G. A., Poelman, D. R., & Dullemond, C. P. 2009, ApJ, 704, 1471, doi: 10.1088/0004-637X/704/2/1471 Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416, doi: 10.1016/j.icarus.2010.12.002 Mu˜noz-Romero, C. E., Oberg, K. I., Banzatti, A., et al. 2024, ApJ, 964, 36, doi: 10.3847/1538-4357/ad20e9 Mulders, G. D., Ciesla, F. J., Min, M., & Pascucci, I. 2015, ApJ, 807, 9, doi: 10.1088/0004-637X/807/1/9 Muller, T. W. A., & Kley, W. 2012, A&A, 539, A18, doi: 10.1051/0004-6361/201118202 Romero-Mirza, C. E. 2023, iris: InfraRed Isothermal Slabs, 0.2.9, doi: 10.5281/zenodo.10369000 Najita, J. R., Carr, J. S., Pontoppidan, K. M., et al. 2013,  \nApJ, 766, 134, doi: 10.1088/0004-637X/766/2/134  \nNeuhauser, R., Guenther, E. W., Wuchterl, G., et al. 2005,  \nA&A, 435, L13, doi: 10.1051/0004-6361:200500104 Picogna, G., & Marzari, F. 2013, A&A, 556, A148, doi: 10.1051/0004-6361/201321860 Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114, doi: 10.1051/0004-6361/201118204 Piso, A.-M. A., Oberg, K. I., Birnstiel, T., & Murray-Clay, R. A. 2015, ApJ, 815, 109, doi: 10.1088/0004-637X/815/2/109 Pontoppidan, K. M., Salyk, C., Banzatti, A., et al. 2024, ApJ, 963, 158, doi: 10.3847/1538-4357/ad20f0 Rafikov, R. R., & De Colle, F. 2006, ApJ, 646, 275, doi: 10.1086/504833 Salyk, C. 2022, csalyk/spectools ir: First release, v1.0.0, Zenodo, doi: 10.5281/zenodo.5818682 Salyk, C., Herczeg, G. J., Brown, J. M., et al. 2013, ApJ, 769, 21, doi: 10.1088/0004-637X/769/1/21 Salyk, C., Lacy, J., Richter, M., et al. 2019, ApJ, 874, 24, doi: 10.3847/1538-4357/ab05c3 Salyk, C., Pontoppidan, K. M., Blake, G. A., et al. 2008, ApJL, 676, L49, doi: 10.1086/586894 Salyk, C., Pontoppidan, K. M., Blake, G. A., Najita, J. R., & Carr, J. S. 2011, ApJ, 731, 130, doi: 10.1088/0004-637X/731/2/130  \nFigure 27. Corner plot for the two-temperature fit to FZ Tau. Dashed lines indicate the 1 σ percentiles. The 2D histogram contours correspond to the 0.5, 1.0, 1.5, and 2.0 σ confidence levels. The histograms have been smoothed by a Gaussian kernel of width 0.02% of the prior span.  \n<!-- image -->  \nSemenov, D., & Wiebe, D. 2011, ApJS, 196, 25,  \ndoi: 10.1088/0067-0049/196/2/25  \nSpeagle, J. S. 2020, MNRAS, 493, 3132,  \ndoi: 10.1093/mnras/staa278  \nToci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, MNRAS, 507, 818, doi: 10.1093/mnras/stab2112 van der Marel, N., & Mulders, G. D. 2021, AJ, 162, 28, doi: 10.3847/1538-3881/ac0255  \nvan Dishoeck, E. F., Bergin, E. A., Lis, D. C., & Lunine, J. I. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 835-858,  \ndoi: 10.2458/azu uapress 9780816531240-ch036 Wells, M., Pel, J. W., Glasse, A., et al. 2015, PASP, 127, 646, doi: 10.1086/682281  \nWoitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383, doi: 10.1051/0004-6361/200911821  \nFigure 28. Same as Figure 27, for the power law column density fit.  \n<!-- image -->  \nWoitke, P., Min, M., Thi, W. F., et al. 2018, A&A, 618, A57, doi: 10.1051/0004-6361/201731460  \nWu, Y.-L., Sheehan, P. D., Males, J. R., et al. 2017, ApJ, 836, 223, doi: 10.3847/1538-4357/aa5b96  \nXie, C., Pascucci, I., Long, F., et al. 2023, arXiv e-prints, arXiv:2310.13205, doi: 10.48550/arXiv.2310.13205  \nZagaria, F., Rosotti, G. P., & Lodato, G. 2021a, MNRAS, 504, 2235, doi: 10.1093/mnras/stab985  \n- -. 2021b, MNRAS, 507, 2531,  \ndoi: 10.1093/mnras/stab2024  \nZhang, K., Pontoppidan, K. M., Salyk, C., & Blake, G. A. 2013, ApJ, 766, 82, doi: 10.1088/0004-637X/766/2/82  \nFigure 29. Same as Figure 27, for the exponentially tapered power law column density fit.  \n<!-- image -->'}