Imprint of urbanization on snow precipitation over the continental USA

Domain: earth_science
Published: 2024-03-15
Editorial outcome: major_revision

Abstract
Urbanization can alter the local climate through modifications in land-atmosphere feedback. However, a continental scale evaluation of its influence on precipitation phase remains unknown. Here, we assess the difference in the likelihood of snow dominated events (SDEs) over 7,415 urban and surrounding non-urban (buffer) regions across the continental United States. Among 4,856 urban-buffer pairs that received at least five SDEs per year, 81% of urban regions are characterized by a smaller snow probability, 99% by a lower frequency of SDEs, and 57% by faster declining trends in SDEs compared to their buffer counterparts. Notably, urban (buffer) regions with lower snow probability are often characterized by higher net incoming and sensible energy fluxes as compared to buffer (urban) regions, thus highlighting the influence of land-energy feedback on precipitation phase. Results highlight a clear imprint of urbanization on precipitation phase and underscore the need to consider these influences while projecting hydro-meteorological risks.

Introduction
Snow precipitation supports water resources for more than 1.5 billion people1. Alarmingly, climate warming is reducing the fraction of precipitation falling as snow2,3. Such a transition is altering water budget partitioning4, flood risks5,6, and water supply7. In addition, changes in snow vs. rain frequency and/or amount have the potential to also impact the timing and magnitude of seasonal and diurnal peak flows8,9, summer baseflows10, evapotranspiration11, nutrient cycling12,13, land–atmosphere feedback14, soil erosion15, vegetation responses16, orographic precipitation generation17, and a host of other coupled ecosystem processes and functions18,19. Especially within the urban regions, the phase of precipitation also influences traffic mobility20 and road accident risks21,22, frequency and intensity of deicing salt application and its consequent impacts on water quality23, contamination from runoff24,25, and vulnerability from extreme hydro-meteorological events26. While the influence of climate variations and changes on the likelihood of snow vs. rain has been studied extensively2,3,27,28, whether, to what extent, and where urbanization impacts precipitation phase over the continental USA (CONUS) remains unknown. Notably, an exclusive focus on evaluating alterations in precipitation phase solely due to global warming, without considering the potential regional impacts of urbanization, would lead to deficient or ineffective projections of its consequences and subsequent mitigation strategies.

Given the wide-ranging consequences of precipitation phase, several of which are noted above, it is imperative to examine the role of urbanization on precipitation phase in order to support long-term water resource planning, ensuring sustainability of ecosystem services, and infrastructure risk assessment and design under changing climate and land cover.

Urbanization’s impacts on surface and air temperature due to the so-called urban heat island (UHI) effect29–32 has been studied extensively, although a majority of these studies have focused on warm/summer periods when snow precipitation is negligible. Furthermore, most past studies on UHI effect did not specifically assess its impact on temperature on the days receiving precipitation, especially with high likelihood of snow. A few studies on the evaluation of the role of urbanization on precipitation phase indicate that urbanization may decrease snow precipitation33–35 because of UHI. Notably, most of these studies have been conducted at local to regional scale. Over the continent, urbanization’s impacts on precipitation phase may however experience significant spatial heterogeneity due to the fact that several locations actually experience urban cooling, instead of heating36,37. In fact, the Global Urban Heat Island DataSet, v1 (2013)38, which consists of data pertaining to differences in averaged daytime maximum and nighttime minimum temperature during summer months of 2013 for 31,500 urban and corresponding buffer regions reveal that 8949 (28%) urban regions during daytime and 11,514 (36%) regions during night are cooler than their buffer counterparts. In addition, given urbanization is also known to alter precipitation temporality and its characteristics39–41, its eventual impact on snow frequency and amount is likely to be influenced by changes in air temperature as well as the temporal distribution and frequency of precipitation events.

This study assesses the differences in modeled average snow probability (SP) between urban areas and their buffer regions, henceforth termed urban–buffer pairs (UBPs), throughout the CONUS. Also evaluated are the differences in the amount of annual precipitation delivered by precipitation events with high snow likelihood. To this end, snow-dominated events (SDEs), i.e., events with likelihood of snow being higher than 50% and 75%, henceforth referred as SDE50 and SDE75, respectively, are identified. We also assess the contrasts in land–atmosphere energy fluxes in urban and buffer regions for each UBP to tease out its potential role on snow precipitation differences between the two regions. Finally, we evaluate the trends in the frequency of SDE50 and SDE75, as well as the snow probabilities associated with them, in all selected UBPs. Unless stated otherwise, reported snow probabilities, their contrasts, and subsequent analysis all throughout this study are based on modeled estimates of likelihood of snow precipitation. More details regarding the selection of UBPs and the data used in this analysis, approach to derive event SP for UBPs, and metrics and thresholds used to define SDEs and assess the differences in their characteristics between urban and buffer regions are presented in “Methods.”

Results
Average SP contrast between urban and buffer regions
Comparison of the average of modeled SP of SDEs between urban and corresponding buffer regions reveals that in about 81% of UBPs (3944 out of 4856), differences in average SP (daSP, see “Methods” for more details) of SDE50s between buffer and urban regions are positive (Fig. 1). In other words, an overwhelming fraction of buffer regions encounter SDEs with higher average SP as compared to their urban counterparts. Relatively high magnitudes of positive and negative daSP are existent mainly in the western part of CONUS, while some also lie in the northeastern US as well. UBPs with negative daSP, though in a small fraction, are spread across the CONUS with greater prevalence in the northern Great Plains region. There are 1387 UBPs (out of 4856) with statistically significant daSP (p<0.05) based on z-score test. A very high percentage (91%) of these UBPs have positive daSP (Fig. S4). Similar dominance of buffer regions with relatively high SP is noted in case of SDE75 as well, where 79% (3578 out of 4553) of UBPs show positive daSP (Fig. S4) and this percentage increases to 84% (851 out of 1013) for UBPs with statistically significant daSP (Fig. S4). Notably, the spatial pattern (areas with relatively high differences and regions with negative daSP) are consistent irrespective of SP thresholds (50% or 75%) for SDEs.

The SP estimates presented above are based on Dai formula (see “Methods: Derivation of event SP for UBPs,” for more details), which is based on air temperature. Evaluation based on Jennings formula, which considers the role of both temperature and relative humidity (RH) on SP, also shows similar results. For example, 86% (84%) of UBPs exhibit positive daSP for SDE50 (SDE75) when calculated using the Jennings’ formula, while this magnitude is 81% (79%) in case of Dai’s formula (see Figs. 1 and S4). Notably, the daSP obtained using the two approaches exhibit strong linear association (Fig. S4b, c). Considering that both the methods yield very similar results, and the Jennings method requires additional data of RH at 4 km × 4 km resolution and hence possibly the uncertainties inherent with it, only Dai’s formula is used for the ensuing analyses.

Factors controlling spatial distribution and magnitude of urban-buffer SP contrast
To understand the spatial distribution of daSP, first the variation of daSP with difference between mean temperatures of buffer and urban regions (Fig. 2a) during 1980–2020 is evaluated. Mean temperatures for urban and buffer regions are respectively obtained by averaging all three-hourly air temperatures during 1980–2020. The correlation between daSP for SDE50 and the difference in mean temperatures between the two regions is found to be r = −0.57. In other words, cooler buffer regions w.r.t. their urban counterparts generally experience higher average SP. Similar analysis carried out with difference in mean critical temperature (see “Methods: Mean critical temperature for more details”), which captures thresholded mean air temperature during cooler periods only, indicates a stronger association (r = −0.68; Fig. S4a). Absolute value of correlation increases from 0.68 to 0.9 when critical temperatures only during the precipitation events (i.e., mean critical event temperature) are used (Fig. S8b). Considering the SP calculated here (Fig. S8b) are based on air temperature (Fig. S3), high correlation evident in Fig. S8b is expected. Nonetheless, larger correlations between daSP and difference between mean critical temperatures (Fig. S8a) and daSP and mean event critical temperatures (Fig. S8b) w.r.t. daSP and mean temperature (Fig. 2a) highlight that only a limited fraction of spatial distribution in daSP can be explained based on the difference in mean temperatures between urban and buffer regions, as its spatial variability is markedly different than that of temperature during the SDEs.

Previous literature32 indicates that difference in air temperature between urban and buffer regions can be due to contrasts in land–atmosphere energy feedback. As shown above, given that the distinctions in daSP between urban and buffer regions are more related to differences of air temperature during SDEs between the two regions rather than overall difference in average temperatures (Figs. 2a and S8), next we assess the differences in land–atmosphere energy exchange only for the days on which SDE50 occurs in either urban or buffer area (hereafter referred to as SDE50-days) of an UBP (see “Data of UBPs and climate” and “Assessment of differences in land–atmosphere energy fluxes between urban and buffer regions” in “Methods” for more details). For each UBP, we first evaluate the percentage of SDE50-days a given energy component is higher in the buffer region compared to the urban. For each energy component, the averages of these percentages over the selected number of UBPs are obtained next. For UBPs where average SP of SDEs is higher in buffer regions (i.e., daSP is positive), SWnet, Rn, and Qh are all greater in buffer regions on an average of 39.78%, 39.67%, and 40.98% of the days, respectively (Fig. 2b and Table S1, Part A). In other words, these energy fluxes are lower in buffer on a majority of SDE50-days. Mean difference (buffer − urban) in the magnitude of these fluxes across all UBPs and SDE50-days also indicate that, on average, these energy fluxes are of lower magnitude (Fig. S9) in buffer regions with higher SP relative to their urban counterparts. In contrast, LWnet is greater in buffer regions on 54.07% of the SDE50-days and Qle is greater in 50.7% of the days. Given the spatial proximity of urban and buffer regions within an UBP, lower SWnet and Rn in buffer regions on a majority of SDE50-days is largely attributable to higher buffer albedo. In fact, 83.2% and 80.83% of SDE50-days with SWnet and Rn being lower in the buffer region also have albedo higher in them. Similar results are noted for the UBPs with negative daSP (Fig. 2b and Table S1, Part B). For example, overall, buffer regions with negative daSP experience higher SWnet, Rn, and Qh on 63.68%, 65.08%, and 66.73% of the SDE50-days, respectively. More than 85% of these days are characterized by higher urban albedo (Table S1, Part B).

Overall, the results reveal that the region (urban/buffer) with relatively high average SP generally experiences lower SWnet, Rn, and Qh, and higher land surface albedo. A higher albedo can decrease net incoming radiation by reducing SWnet. While LWnet increases due to concomitant decrease in surface temperature which in turn reduces the outgoing longwave radiation, but given that change in LWnet is generally smaller than SWnet, Rn still is lower with higher albedo. The decrease in Rn and surface temperature consequently diminishes sensible heat feedback42 to the atmosphere. However, the extent of impact is also determined by several other mediating factors32, including land roughness properties, emissivity, moisture sources, etc. Since Qle approximately equals Rn minus Qh, and as both Rn and Qh are lower in regions with higher SP, Qle shows no marked contrasts between the two regions (Figs. 2b and S9), with it being higher on only around 51% of the SDE50-days across the UBPs. Overall, considering that Qh has a strong influence on air temperature43, the contrasts in it between urban and its buffer counterpart contribute to differences in SP. In fact, among the UBPs with higher average SP in buffer regions, i.e., UBPs with positive daSP, 78.64%, 78.27%, and 78.11% of the days with lower Qh in buffer have lower Rn and SWnet, and higher albedo, respectively (Table S2, Part A). Corresponding fractions of 80.28%, 79.59%, and 85.01% for UBPs with higher average SP for urban regions further attest the findings (Table S2, Part B).

It remains unclear whether the increased albedo that leads to reduced Rn and Qh, subsequently resulting in lower air temperatures and an increased likelihood of snowfall, is primarily due to differences in snow cover conditions or land cover characteristics between the two regions. To assess this, we next evaluate the energy fluxes for SDE50-days but only when both urban and their buffer counterparts have snow-free ground (see “Methods” for the details about snow-free ground conditions). For UBPs with positive and statistically significant daSP (i.e., in about 81% of UBPs with statistically significant daSP), albedo in urban is lower for snow cover-free days (Table S1, Part C). However, for UBPs with negative daSP, the fraction of days when buffer Rn and Qh fluxes are lower than urban counterparts are reduced for snow cover-free conditions (Table S1, Part D). Similar evaluations of energy contrasts for snow-covered days (Fig. S10 and Table S3) further reconfirm the underlying role of albedo differences between urban and buffer regions on SP differences. These results highlight that albedo disparities, both during snowcover-free and snow-covered days, are closely associated with disparities in Rn and Qh between urban and buffer counterparts. The overall association between positive daSP and disparities in energy fluxes between urban and rural regions should not be construed as implying that every SDE occurrence will experience these associations. This is particularly true as other dynamic factors such as land surface properties and moisture sources, may also play a mediating role in determining the discrepancies of energy fluxes between urban and rural regions. Another reason is that energy fluxes at daily scale can only be considered as an approximate descriptor of air temperature during SDEs, especially if a day receives multiple SDEs. To assess the imprint of buffer land covers on daSP contrasts on SDE50 days with snow cover-free ground conditions, we stratify the aerial fraction of 11 land covers: water, developed areas, deciduous forest, evergreen forest, mixed forest, shrub (open), grassland, hay, cultivated crops, woody wetlands, and herbaceous wetlands. Fractions of each land cover are first obtained by mapping National Land Cover Database land cover to each buffer region with statistically significant daSP (Fig. S4). The average of these fractions for each land cover is then evaluated separately over UBPs with positive and negative daSP. Results (Table S4) indicate that UBPs with negative daSP are distinguished by a higher forest coverage (deciduous, mixed, and evergreen forests). In contrast, UBPs with positive daSP have a higher proportion of cultivated crops. This imprint is also reflected in sparse occurrence of negative daSP UBPs in croplands of mid-western US. While forest (cultivated crop) land covers usually are characterized by lower (higher) albedos than in built-up areas, it is to be noted that land–atmosphere energy feedback in buffers are also determined by differences in surface and aerodynamic conductances44 as well. Moreover, the impacts of additional land covers, though individually minor, could wield a more substantial influence than the predominant land cover in the area. Mapping the influences of each land cover on different energy components with desired certainty currently remains beyond the scope of this study.

Additional analysis is carried out to assess the characteristics of SDEs that result in positive daSP over the continent. To this end, first the SDE50 falling within specific SP ranges are identified in both urban and buffer regions. Next, the difference in average SP for each of these ranges is evaluated for each UBP. Average of these differences (Fig. 2c) overall UBPs for different SP ranges indicates that the largest difference in SP is for events falling within the 90–100% SP range. In other words, buffer regions generally have much higher SP than urban regions for SDEs with SP of 90–100%. For other SP ranges, the differences are relatively much smaller. To further confirm that SDEs with SP in the range 90–100% are the main driver of positive daSP overall, we next assess whether these events have at least similar likelihood of occurrences as other SP ranges. To this end, we evaluate the ratio of SDEs with SP falling within a specified range to the total number of SDE50 in each urban or buffer region. Next, the difference of these ratios between buffer and urban regions of each UBP is calculated. Average of these differences overall UBPs (Fig. 2d) show that for the SP interval 90–100%, the fractions of SDEs between buffer and urban have the higher difference. These results illustrate the predominant influence of SDEs with SP of 90–100% on why the average SP is generally higher for buffer region. Results for SDE75 agree with those for SDE50 (Fig. S11).

Is there an urban–buffer contrast in precipitation amount delivered by SDEs?
Although an overwhelming fraction of UBPs reveal higher SP in buffer regions, only 54% of the buffer regions, out of 4856 UBPs, receive a higher amount of annual average precipitation delivered by SDE50 (Fig. 3a). This percentage reduces to 52% when urban and buffer regions are compared on the basis of annual average percentage of precipitation received through SDE50 (Fig. 3b). Notably, 99% of buffer regions experience a higher frequency of annual average SDEs (SDE50) (Fig. 3c). Despite this, the reason an overwhelming number of buffer regions do not encounter higher amount of precipitation (or even higher fraction of snow to total precipitation) than their urban counterparts is because for 97% of the UBPs, the buffer regions receive SDE50 at lower intensities as compared to their urban regions (Fig. 3d). In other words, more number but smaller intensity SDEs result in balancing out of the total precipitation delivered by SDEs between urban and buffer regions. Notably, a smaller intensity of precipitation events, though not for SDEs, in buffer regions w.r.t. urban areas has also been reported previously45. Evaluations for SDE75 also show similar results (Fig. S12).

Differences in temporal trends of SDEs between urban and buffer regions and their controls
The trends of annual frequency of SDEs and annual average SP of SDEs are evaluated next (Fig. 4). In case of annual frequency of SDE50 (Fig. 4a), ~74% (3578 out of 4856) UBPs have statistically significant (s.s.) trends (p<0.05) in both urban and buffer regions. ~98% of these UBPs with s.s. trends (3500 out of 3578) exhibit a declining behavior, i.e., the frequency of annual SDE50 decreases over time in these UBPs. Notably, in ~57% of UBPs with s.s. negative trend (1991 out of 3500), urban regions show a more pronounced negative trend in annual frequency of SDE50 compared to buffer. Trend analysis for SDE75 reveals similar findings (Fig. S13a). While these UBPs are distributed all across the CONUS, those with largest discrepancy in trends between urban and buffers either lie in the northernmost states east of the Mississippi River or are nestled in the mountainous regions (e.g., the Rockies of UT and CO, or the Sierras of California) of the western US (Fig. S14). The more pronounced negative trend in SDEs frequency in urban areas, compared to buffers, may result from two factors: (1) a higher temperature increase in urban regions, which can potentially transition SDEs into rain events, reducing their frequency, and (2) dwindled precipitation events in urban areas, which include a reduction in SDE frequency as well. To validate these possibilities and identify the dominating factor affecting the evolution of annual frequency of SDEs, trends of mean annual temperature and annual frequency of precipitation events are obtained for the UBPs with urban regions showing pronounced negative trends (Table S5, Parts A and B). Among the 1991 UBPs where annual frequency of SDE50 exhibits a more negative trend than buffer, a noteworthy 96% of them (1908 out of 1991) show statistically significant negative trends for annual frequency of precipitation in both urban and buffer regions, with the urban regions showing a more substantial negative trend (Table S5, Part A). In contrast, only 40% of the 1991 UBPs have higher magnitude of annual average temperature trends for urban regions. A similar pattern is noticed while examining trends for SDE75 (Table S5, Part B). These results indicate that for urban regions with higher negative trend of SDE frequency, the primary factor for the expressed contrasts is the disparity in the trends of annual precipitation frequency, rather than the air temperature.

For annual average SP, a smaller fraction (~20%) of UBPs show s.s. trends (982 out of 4856) in both urban and buffer regions. However, 98% of these UBPs (958 out of 982) exhibit declining trends in both urban and buffer regions. In 70% of these UBPs (675 out of 958), urban regions show a more pronounced negative trend than buffer areas. Trend analysis for SDE75 reveals similar findings (Fig. S13b). Largest differences in trend magnitude are noted in central and northwestern US (Fig. S14c, d). Next, we assess the potential impacts of differing trends in air temperature and critical event temperature, the two variables that show a significant association with the difference in average SP between urban and buffer regions (Figs. 2a and S8), on the varied trends in average SP. The analysis reveals that, more so than the trend in annual average temperature, higher positive trends in annual average event critical temperature are closely associated with a more rapid decline in the annual average SP of SDEs (Table S5, Parts C and D).

Methods
Subsequent sections provide comprehensive information pertaining to the data and the adopted procedure for evaluating SP contrasts between urban and corresponding buffer regions (Fig. S1). This includes details regarding delineating buffer regions around each urban area, downscaling and spatially explicit mapping of climate variables, obtaining spatially-averaged representative temperature and precipitation time series for both urban and buffer regions, and finally obtaining differences between average snow probabilities for SDEs (SP > 50% and 75%). Additionally, these sections elaborate on the metrics and thresholds employed for the evaluation of the differences in snow precipitation characteristics between urban and buffer regions.

Data of UBPs and climate
In this study, the Global Rural-Urban Mapping Project, Version 1 (GRUMPv1): Urban Extent Polygons, Revision 02 is used to obtain the data of urban regions38. GRUMPv1 provides data of 75,445 urban areas across entire globe out of which 9359 urban regions are situated within the CONUS (Fig. S1a). The data have been widely used in studies to map and understand UHIs50–52. Corresponding to each urban area, surrounding non-urban or buffer zone is mapped using the Buffer function in ArcGIS Desktop 10.6.1 (Fig. S1b). Buffer corresponding to each of the urban area is drawn using a 10 km distance threshold from the boundary of demarcated urban areas. It is to be noted that while different approaches exist for the selection of size of buffer regions, including those based on distance thresholds (e.g., use of 1 km threshold37 and 10 km threshold32), fraction of urban area-based approach53,54 (e.g., area of buffer region between 50% and 150% of urban area), and pixel-based approach55, the choice of distance threshold used here is guided by the need to feasibly and representatively ascertain SP contrasts within UBPs over the entire CONUS. Use of 10 km threshold allow sample samples from the buffer region. Increasing the width of buffer region beyond 10 km engulfs surrounding urban regions while reducing the width to say 1 km is not advisable as the downscaled climate variables are at 4 km resolution. Notably, the 10 km threshold has also been selected for analysis in other previous studies32,38 related to UHI effect as well.

For ensuing analysis, precipitation and air temperature data from the NLDAS project56 are used. These data are available at 0.125° spatial and an hourly temporal resolution, and have been widely used in a range of disciplines such as hydrology57,58, UHI impact on public health59,60 and agriculture61. Although other data at much finer spatial resolution exists, they are however generally at a temporally coarser resolution. The sub-daily resolution of the temperature and precipitation time series in the NLDAS dataset is crucial for appropriate apportioning of precipitation into snow and rain62 and also for evaluation of SP46 for a given event.

To evaluate the energy exchange between land and atmosphere on the days when SDE50s are encountered by either of the two regions in an UBP (SDE50-days), outputs from the Noah land surface model simulations, which are also obtained for the NLDAS project56, are used. The Noah model serves as the land component in several coupled land–atmosphere modeling systems, including the Weather Research and Forecasting (WRF) regional atmospheric model, the NOAA-NCEP coupled Climate Forecast System (CFS), and the Global Forecast System. The energy components that are considered for this analysis include (1) net shortwave radiation flux (SWnet), (2) net longwave radiation flux (LWnet), (3) net radiation flux (Rn), (4) sensible heat flux (Qh), and (5) latent heat flux (Qle). Net radiation flux is calculated as the sum of net shortwave and net longwave radiation fluxes from the Noah model output. Data of surface albedo from the model output are also used.

Derivation of event SP for UBPs
The procedure to estimate the probability of precipitation event being snow involves multiple steps (Fig. S1). These include:

Climate data mapping onto 4 km × 4 km grid over the CONUS. Considering that more than 80% of the 9359 urban regions situated within the CONUS, based on GRUMPv1 data, have areas less than the spatial extent of a NLDAS grid (area ~0.125° × 0.125° or ~190 km2 at equator), the coarse resolution temperature data are downscaled to 4-km grid using inverse distance bilinear interpolation approach. For obtaining interpolated temperature at any given 4 km × 4 km grid center, the distances between the center of this grid and centers of four surrounding NLDAS grids are used as weights. Similar spatial downscaling of temperature data has also been performed in other data products63.

Temporally, temperature data are up-scaled to 3-h resolution. This is because the parameters of SP function used in this study (explained later in the “Methods: SP in urban and buffer regions” section) were established using data at 3-hourly resolution. Precipitation data from NLDAS are, however, transferred as-is to the encompassing 4 km × 4 km resolution grids, in part because precipitation generally shows significant spatial heterogeneity governed by additional meteorological, physiographic, and topographic attributes64–66. Adopting a simple interpolation scheme for precipitation data may introduce large uncertainties leading to erroneous findings56. Notably, also within the reanalysis datasets67, spatial distribution of precipitation are generally determined by the model outputs or spatially-explicit radar based observations, instead of using simple (e.g., bilinear) spatial interpolations of the observations.

Obtaining spatially averaged climate of urban and buffer areas. For each UBP, all 4 km × 4 km grids with more than 50% of their areas within the urban or buffer regions, respectively, are identified (see Fig. S1c). Out of 9359 UBPs (Fig. S2a) within the CONUS, 7415 UBPs (Fig. S2b) have at least one 4 km × 4 km grid each within urban and buffer regions. For the ensuing analyses, these 7415 UBPs are selected. For computational tractability, a single temperature and precipitation time series are obtained by spatially averaging data of selected grids in each of the 7415 urban and buffer regions. For example, an urban region consisting of n number of grids (at 4 km resolution) with more than 50% of their areas within the region, temperature time series for the region is generated by averaging the temperature time series for n grids.

SP in urban and buffer regions. Phase of precipitation at the ground surface is influenced by the properties of the atmosphere through which it travels. Different models for precipitation phase determination68 have been proposed in previous studies, including those based on (1) air temperature and RH thresholds46,69,70 and (2) atmospheric models with microphysics schemes that track a hydrometeor from its formation71. Here we obtain SP based on a hyperbolic tangent function (Fig. S3) whose parameters have been established using 3-hourly weather reports data46. The function (Eq. 1) has been used in numerous studies to perform snow vs. rain partitioning27,47. Hereinafter, this function will be referred as Dai formula.

SP = 48.2292 × (tanh(0.7205 × (temperature − 1.1662)) + 1.0223) (1)

SP is also evaluated using a bivariate logistic regression model47 (termed as Jennings formula henceforth) that links SP with both RH and temperature

SP = 100 / (1 + exp(α + βTs + γRH)) (2)

where α = −10.04, β = 1.41, γ = 0.09. For this, the RH is calculated based on Clausius–Clapeyron equation at the spatial and temporal resolutions of NLDAS data using Eq. 3.

RH = 0.263 × specific humidity × pressure / exp((17.67 × (Ts − 273.16)) / (Ts − 29.65)) (3)

where specific humidity is expressed in kg/kg, pressure in Pascal, and Ts in Kelvin. Specific humidity and pressure data are obtained from NLDAS datasets. RH obtained in previous step is spatially downscaled (at 4 km resolution) and temporally up-scaled at 3-h resolution using the same approach as that for the temperature (see “Climate data mapping onto 4 km × 4 km grid over the CONUS” in “Methods”).

Metrics and thresholds to assess the differences in snow precipitation characteristics between urban and buffer regions
As the premise of this study is to evaluate the impacts of urbanization on snow events, here we focus on SDEs, i.e., events for which the SP ≥ 50%, i.e., the probability of event being in solid phase is higher or equal to than it being in the liquid phase. We also define the metrics and variables that allow the assessment of contrasts in snow characteristics between urban and buffer regions. Finally, UBPs where robust statistical analyses of the variables under consideration can be performed are identified.

Thresholds used to define SDEs. Here two threshold snow probabilities are used to define SDEs. These include SDE with SP greater than 50% and 75%. These events are termed as SDE50 and SDE75, respectively. These snow probabilities correspond to air temperature (Tthresh) = 1.146 °C and 0.341 °C, respectively (Fig. S3), based on the Dai formula.

Metrics used to assess snow contrasts between urban and buffer regions. Six metrics are used to assess the impact of urbanization on snow characteristics. The first metric used here is the difference (buffer–urban) in average snow probabilities for SDEs, henceforth called the daSP, between the buffer and the urban regions during 1980–2020 (Fig. S1d). This metric helps identify whether the buffer or the urban region receive SDEs with a higher average probability of snow. The second metric quantifies the difference (buffer–urban) in the annual average amount of precipitation delivered by SDEs during the analysis period, i.e., 1980–2020. The third metric evaluates the difference (buffer − urban) between the annual average percentages of precipitation that is delivered by SDEs. The fourth metric compares the intensities of precipitation (=precipitation amount/number of SDEs) delivered by SDEs. It is quantified in mm/events/year by subtracting the average annual intensity of SDEs over the urban region from that in the buffer. The final two metrics compare the temporal trends of (1) annual frequency of SDEs and (2) annual average SP of SDEs (for both SDE50 and SDE75) for each UBP. The trends are calculated using the non-parametric Sen Slope estimator.

Mean critical temperature. Mean critical temperature and mean critical event temperature are evaluated to assess the role of air temperature during cooler periods on SP. Mean critical temperature for urban/buffer region is obtained by only averaging the air temperature values below 1.146 °C during 1980–2020, after substituting the temperature data below −4 °C with −4 °C. The reason for using 1.146 °C as an upper threshold is that this temperature corresponds to 50% SP (Fig. S3). As SP below −4 °C remains almost constant, the temperature values below this threshold are replaced by −4 °C. The selection of these thresholds enables evaluation of a temperature metric, which can be easily calculated and is expected to capture the variations of SP across regions. Mean critical event temperature is calculated in a similar manner as the mean critical temperature. However, temperature values only during the precipitation events are considered for its evaluation.

Selection of UBPs. Even though the SPs are obtained for all 7415 UBPs, all evaluations in this study are performed only for the UBPs that receive at least 5 SDEs per year in both urban and buffer regions. Out of 7415 UBPs (see Fig. S2b), 4856 (4553) UBPs encounter at least 5 SDE50 (SDE75) per year (Fig. S2c, d) in both urban and buffer regions of a UBP. This additional criterion is imposed to ensure sufficient sample size (frequency of SDEs) while obtaining the metrics.

Assessment of differences in land–atmosphere energy fluxes between urban and buffer regions
Spatially averaged time series of energy components and albedo for each UBP are obtained using the same procedure as temperature as discussed above in the sections “Climate data mapping onto 4 km × 4 km grid over the CONUS” and ”Obtaining spatially averaged climate of urban and buffer areas.” These data are obtained for a 41-year period (1980–2020).

Contrasts in energy fluxes are evaluated for 1387 UBPs with statistically significant daSP (based on 41 years of SP calculations for SDE50s (see supplementary Fig. S4)). For each of these UBPs, spatially and temporally averaged daytime (i.e., hours when net shortwave radiation flux > 0) energy fluxes such as SWnet, LWnet, Rn, Qh, Qle, and albedo are calculated individually for both urban and buffer regions. The evaluations are performed for all days receiving SDE50 in either urban or buffer region. As antecedent states at sunrise and cumulative Qh to the atmosphere drives the net change in boundary layer thickness and air temperature43 during a day, daily energy fluxes between urban and buffers are compared. Evaluations are also performed separately for a subset of these days when (1) the ground is snow cover-free and (2) ground is covered with snow, for both the urban and buffer regions. Snow cover-free days based on NLDAS data are defined as the days when modeled albedo of both the regions is <0.35 and snow-covered days are identified by complimentary magnitude of albedo (≥0.35). The threshold is conservative given that the upper limit of all land covers considered in the NLDAS project has an albedo of 0.33.

Data availability
NLDAS gridded precipitation, temperature, energy fluxes, and albedo hourly data can be downloaded from https://disc.gsfc.nasa.gov/datasets?keywords=NLDAS. The Global Rural-Urban Mapping Project, Version 1 (GRUMPv1): Urban Extent Polygons, Revision 02 can be retrieved from https://sedac.ciesin.columbia.edu/data/set/grump-v1-urban-ext-polygons-rev02/data-download. Data used for the analysis of snow characteristics contrast between urban and buffer regions may be accessed from https://zenodo.org/records/10723791. Source data are provided with this paper.

Code availability
The codes to generate figures may be accessed from https://zenodo.org/records/10723791. Detailed information about required input files and steps to run the codes are provided in readme.docx.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript explored the differences in the snow dominated events(SDEs) between urban and\nsurrounding buffer regions across the confinental United States. The main conclusion was that the\nurbanizafion has a negafive impact on the occurrence probability of snow events, indicafing that urban\nregions tend to have smaller snow probability, lower frequency, and decline in SDEs. These findings are\ninteresfing. However, the following issues should be clarified.\nScienfific relevance:\nThe major issue is that the current manuscript does not provide a very good scienfific quesfion, and it\nsuffers from limited physical explanafions, and an overall lack of discussions as how this study builds\nupon the large body of studies around the world.\n(1)Do the authors consider the changes of urban regions during urbanizafion process?\n(2)Could the authors explain what does posifive daSP(in line 91) and other negafive mean? Could please\ngive a clear explanafion?\n(3)Is it reasonable to using a constant 10km distance threshold for all the urban regions?\n(4)The methodology, especially the detailed informafion of the climate data used should be furthered\nclarified.\n(5)As the authors menfioned, there are sfill some uncertainfies in the manuscript. Therefore, how did\nthe authors do their best to minimum?\n(6)The first sentence in line 216 was too absolute. Did the authors review all the previous literature?\n(7)It is strongly recommended to proofread the paper thoroughly again.", "reviewer_2": "The manuscript documents differences between urban and adjacent buffer zones in snow to\nprecipitafion rafio, snow frequency, snow intensity, and trends in these metrics. The noteworthy results\ninclude that a) buffer areas tend to receive more frequent but less intense snowfall events that are\ncharacterized by a higher snowfall fracfion relafive to urban areas and b) trends in the frequency of snow\nevents and their snow to rain rafio are more negafive in urban areas than buffer areas.\nThe work is original and the claims are supported by the evidence. However, the broader impacts of the\nwork and its significance are relafively weak, at least as currently framed in the paper. The authors\nmenfion why precipitafion phase mafters, but stronger mofivafion for looking at effects of urbanizafion\nspecifically on precipitafion phase would improve the impact of the paper. While I believe the work is\nabsolutely worth publishing, I am not sure the significance is great enough for publicafion in Nature\nCommunicafions.\nI had several quesfions about the methods, mostly clarificafions of text that I found unclear. My only\nmajor concern with the methods, based on my current understanding, is the spafial resolufion and\ndownscaling of the climate data. I am concerned that inverse distance weighfing downscaling of ~12km\ndata to 4km does not capture the real spafial heterogeneity of climate fields at 4km, rather it just\nsmooths a coarse grid. I was concerned that this would make it difficult to capture contrasts between\nurban-buffer pairs, but the analysis demonstrated contrasts, so perhaps this is less of a problem that I\ninifially thought. Addressing this would require using a different climate dataset (and I understand 4km\nhourly data over CONUS is a big ask) or a more sophisficated downscaling method. The lafter might be\nworth considering.\nOverall, the paper is thorough, interesfing, and mostly well presented.\nSincerely,\nAbby Lute\nA few less important but general comments:\nI’m not super familiar with Nature Communicafions formafting. Is it normal to have such lengthy pointers\nto the methods secfion in the main text and to have so much content in the figure capfions?\nIt would be very helpful (and more pleasant to look at) if the color coding of posifive and negafive daSP\nin the figures and tables throughout the paper was consistent.\nI highlighted some grammafical errors, but there were many that I didn’t highlight. Feedback from a\nscience writer or editor would improve the readability of the paper.\nBelow I provide more specific feedback.\nAbstract\nLine 16- ‘the’ not needed\nLine 16- The fact that it is unknown is a relafively weak jusfificafion for a study. Why do we want to know\nit? Why does it mafter if precipitafion phase depends on urbanizafion? Answering these quesfions might\nprovide a stronger mofivafion for the study and get the reader more interested.\nI don’t think you need so many numbers in the abstract, it makes it a liftle harder to get at the big\npicture. Especially in Line 19, these numbers could be omifted for readability.\nLine 19- delete ‘an’\nLine 20- could ‘relafively rapid decline’ be rephrased to clarify that we are talking about long-term\nhistorical trends here?\nLine 22- Suggesfion: replace ‘are found to be usually’ with ‘are often’\nLines 21-24- is this lower of snow probability in general or relafive to their counterpart (urban/buffer)? It\nis unclear at the moment. And it seems like ‘fracfion’ could be left out.\nLines 24-26- Rephrase. Precipitafion phase is not suscepfible to urbanizafion, it is influenced by it\nthough. It might be fairer to say ‘Results show that trends in snow frequency and snow to rain rafio are\nmore negafive in urban areas than buffer areas, highlighfing the effects of urbanizafion on precipitafion\nphase.’\nMain\nLine 29- resource should be plural\nLine 32- have not has\nLine 37- suggesfion: replace ‘does the urbanizafion impact’ with ‘urbanizafion impacts’\nLines 35-38- again this is very weak mofivafion for the study. Why do we care about the effect of\nurbanizafion specifically on precipitafion phase?\nLine 42- replace ‘may result in decrease in’ with ‘may decrease’\nLines 50-51- This sentence is unclear and seems overly wordy. Is the intenfion: “The extent/intensity of\nthe UHI effect has previously been mapped for summer and winter days.”?\nLine 52- 54- check grammar in this sentence\nLine 59- delete ‘as’\nLines 65-66 Just saying trends would be clearer than temporal evolufion.\nAverage snow probability contrasts between urban and buffer regions\nFigure 1- Do I understand correctly that a daSP of 0.5 for example means that there is less than 1%\ndifference in snow probability between urban and buffer?\nFigure S7- please add 1:1 line on panel b) and c). Even though it is mostly obvious, would be nice to add\na panel like a) but for the Dai SP, just to make it easy to compare the two methods.\nFactors controlling spafial distribufion and magnitude of urban-buffer snow probability contrast\nLine 102- is ‘mean temperature’ the mean annual temperature? Mean temperature during SDEs? Please\nprovide a liftle more descripfion. It is hard to interpret the findings presented in this paragraph without\nhaving a liftle befter understanding of the metric.\nFigure S8- average event crifical temperature should be defined in the methods secfion. The gridlines are\ninconsistent in different areas of the plot, which seems odd.\nLine 111- replace Figure references with metric references. E.g. “larger correlafions with average crifical\ntemperature and average event crifical temperature compared to mean temperature highlight…’\nOtherwise the reader has to translate what the figures say in order to make sense of the sentence.\nFigure 2b- This is prefty challenging to interpret as there is a lot of translafion that has to occur. It’s hard\nto suggest alternafives without seeing the data, but I wonder if a barchart showing the mean differences\nin each flux (W/m2) on SDE50 days would work. For each flux there would be two bars- one for posifive\ndaSP and one for negafive daSP. I think it would help to keep the calculafion the same between posifive\nand negafive daSP (i.e. either urban – buffer or buffer – urban, not switching like you’re doing now). This\nwould have the added advantage of showing the relafive magnitudes of the differences in the different\nfluxes. Also, Qle is not defined in the text.\nQle presents different behaviour than the other fluxes (Figure 2b) but this divergent result is not\ndiscussed in the text. How is this behaviour explained?\nLine 151- I don’t think these other factors were noted above.\nLine 157- Why not use energy fluxes during the events instead of for the full day? Maybe this was in the\nmethods and I missed it.\nFigures 2c and d- please remove grey space so that plot content can be enlarged and easier to see. Might\nbe helpful to label the colorbars to indicate what metric they are showing. I prefer a color bar that\ndiverges at 0 for a metric like this. As it is, the colors suggest that green and brown cells mean different\nthings, when in reality they might both be showing posifive values. Why are some heat map cells divided\nby a black line and others aren’t? Same comments apply to Figure S9.\nLines 160-161- by ‘difference in SP’ do you mean the same daSP that’s been discussed throughout the\npaper? If not, how is this different? What is the SP for a difference in SP? For example, if this is daSP then\nis the SP associated with it the urban SP or the buffer SP?\nLine 162- what is the difference in fracfions of SDEs? What is the fracfion that is referred to? This whole\nparagraph should be befter explained in the methods. For example, the grouping into seven bins should\nbe in the methods.\nLine 165- what are frequency-fracfions? Where was this defined?\nLines 159-171- I read this paragraph many fimes and could not figure out what it said. There seem to be\nseveral metric definifions missing, or maybe they are defined in the methods but the same terminology\nis not used here. I was unable to evaluate the science in this paragraph because I couldn’t understand\nthe presentafion of it.\nIs there an urban-buffer contrast in precipitafion amount delivered by snow-dominated events as well?\nLine 174- delete ‘as well’\nLine 175- I would rephrase the first few sentences of this paragraph to say that there is liftle contrast in\na) the amount of precipitafion delivered or b) the precipitafion fracfion between urban and buffer,\ninstead of focusing on the relafively small differences you found.\nFigure 3b- if this is a difference between two fracfions that each range from 0 to 1, how can there be\nvalues greater than 1 or less than -1?\nFigure 3- again the grid lines are inconsistent. Not a huge deal, just annoying. In the capfion, SD50 should\nbe changed to SDE50 and either iSDE needs to be changed to match the abbreviafion in the methods\nsecfion (ipptSDEs) or vice versa.\nLine 187- was there any explanafion for the difference in intensity offered by this previous work?\nDifferences in temporal trends of snow-dominated events between urban and buffer regions, and their\ncontrols\nLine 192- delete both instances of ‘or’\nLines 195-197- This sentence is mostly talking about afSDE and Figure 4a, but then it says ‘both afSDE\nand aaSPSDE decrease with fime in these UBPs’. This sentence makes sense with Figure 4a, but\ncontradicts Figure 4b.\nFigure S11- The clustering in these plots seems very strange. Is this because only stafisfically significant\ntrends are shown? If so, the fact that only these trends are shown should be noted in the capfion. Even if\nthat is the case, I don’t understand why there would be clusters at the 0 temperature trend lines. Trend\nunits should be noted on x and y axis labels or in capfion.\nFigure 4- the descripfion for the smallest semicircle in a) and b) is unclear. What does ‘relafively rapid’\nmean? How about something like ‘trend in urban region is more negafive than buffer region’ or ‘urban\ntrend is more negafive than buffer trend’? Same applies to Figure S12.\nCan you comment on the spafial distribufion of the s.s. trends and of the locafions where urban trends\nare greater than buffer trends?\nConclusions\nLine 216- should be ‘analysis’ not ‘analyses’\nLine 218- rephrase as ‘urbanizafion decreases snow probability’. ‘negafive impact’ sounds like a value\njudgement\nLine 218- ‘Alarmingly’ seems inappropriate here without being explicit about why this should be\nalarming. Also, why would this be? Is it because of interacfions between urbanizafion and changing\nclimate, because characterisfics of urban areas are changing, or something else?\nLines 224-239- I understand the issues of scope and data limitafions, but it seems like one could\ncomment more specifically on potenfial contribufions to urban-buffer differences. For example, it seems\nlikely that buffer regions have more vegetafion and soil surfaces as opposed to concrete/asphalt surfaces\nand this could affect local convecfion, moisture availability, and latent heat flux in specific ways.\nWhat are the implicafions of less snow in urban areas? In many ways this seems like a good thing- it\nwould likely contribute to less congesfion, fewer accidents, and less funding needed to be spent on road\nclearing and repairs. On the other hand, if this is precipitafion is coming as rain instead perhaps it\nincreases flood risk and poses challenges to drainage systems in urban areas which generally have\nsurfaces with lower permeability.\nWhy do buffers tend to get more frequent but less intense SDE’s? What are potenfial explanafions for\nthis? This seems like one of the more interesfing findings to me.\nHow do you explain large negafive daSP values in the West? Is there something unique about these\nlocafions?\nMaybe also menfion that this study does not assess snow on the ground or how snowfall persists on the\nlandscape, but is focused on snowfall. A nice future work piece might be looking at complementary data\nsources (ground observafions or radar) to see both if observafions backup what you see with\nmodeling/reanalysis and also to consider contrasts in snow on the ground / snow persistence.\nMethods\nData of UBPs and climate\nFigure S1- why are there extra curvilinear lines in the rasterized buffer example? Urban legend label\nneeds to be adjusted to be next to the legend key. The map in a) could use a country outline since the\ncolor for the country and the plot background are the same. In b) it looks like there are some excepfions\nto the 10km buffer, where areas within 10km are excluded. Please explain these in the methods.\nA quick glance suggests that the data used for the energy fluxes also contains hourly temperature and\nprecipitafion. Why are these not used instead of the NLDAS data? Also, why use shortwave and net\ninstead of shortwave and longwave?\nDerivafion of event snow probability (SP) for UBPs\nDownscaling of temperature (and energy fluxes) to 4km using an interpolafion approach seems like it\nwould not capture urban/buffer contrasts well at all, especially since more than 80% of urban areas\nconsidered have areas smaller than the NLDAS grid.\nLine 297- the citafion referenced here used inverse distance weighfing to get a grid from stafion\nobservafions to compare to a downscaled temperature dataset, the approach was not used to\ndownscale temperature per se. Perhaps there is a befter citafion you can use to jusfify your approach.\nHow is an event defined? Is it one three-hour fime period with snow? This is important to note since in\nother works snow events have been defined as three days, or other lengths of fiem.\nMetrics and thresholds to assess the differences in snow precipitafion characterisfics between urban and\nbuffer regions\nLine 346- could rephrase as ‘events for which the SP is >= 0.5’\nLines 352-353- Do you mean the snow probabilifies correspond to these temperatures using the Dai\nmethod? Because for the Jennings approach SP would be a funcfion of RH as well and you wouldn’t be\nable to define a single temperature for each SP. Please clarify in the text. Also the capfion of Figure S3\nshould note that this is based on the Dai formula.\nLine 354- please include units for each metric\nLine 356- please add (buffer – urban) or (urban – buffer) to clarify how daSP is calculated.\nLine 360- for pptSDE is this the average precipitafion of an event or the annual sum of precipitafion from\nthese events?\nLine 362- how is the intensity of SDE defined?\nLine 364-367- why are the trends evaluated differently? Why not use the same method for both metrics?\nFor the trend metrics are the annual frequency and annual average SP summed across the urban and\nbuffer components of each UBP? Or averaged?\nThe descripfion of the average crifical temperature makes no sense to me. Why are temperature values\naveraged if they are below 1°C? Are these temperature values at any point in the year or just during\nSDEs? Why do temperature values below -4°C need to be reset to -4°C? How is this a crifical\ntemperature, what does it tell you? Maybe I’m missing something here because this is very unclear to\nme.\nLine 374- is this 5 between the urban and buffer areas, or 5 each?\nAssessment of differences in land-atmosphere energy fluxes between urban and buffer regions\nFigure S4- Suggesfion: I prefer when the color ranges are symmetric around 0. I think the text ‘UBPs with\nposifive daSP: ~91%’ belongs in the manuscript text, not embedded in the figure. Histogram needs a y\naxis label. I would also add state outlines below the points to help the viewer befter understand the\nspafial distribufion of the points. And you might remove the box around the histogram plot. All the same\ncomments apply to Figure 1, Figure S5, and Figure S6 as well.\nLine 386- capitalize figure\nLine 389- so this analysis uses days of SDEs, not strictly the 3-hour fime periods when they occur? Why\nnot be more specific about the fime?\nLine 391- replace ‘correspond to’ with ‘are defined as days’ to clarify that this isn’t fact, it is how you’re\ndefining them for the purpose of this study.", "reviewer_3": "The analysis presented uses temperature and precipitafion data from the North American Land Data\nAssimilafion System, combined with a simple stafisfical model of precipitafion type, to evaluate\ndifferences in the probability and quanfifies of snowfall in Urban versus Buffer (adjacent non-urban)\nareas (termed Urban-Buffer pairs). The authors find that Urban areas, owing to lower land surface\nalbedo and associated greater quanfifies of net solar radiafion, net radiafion, and greater sensible heat\nflux, have lower snowfall probabilifies than adjacent buffer areas. The analysis is interesfing, well\nconsidered, and thorough and the contribufion will be valued by many in the scienfific community and\nby the general public. My comments are intended to address a few potenfial issues inherent in the\nanalysis and text.\n1) The authors rely on a model and not observafions. I assume there is a good jusfificafion for this as it\nseems obvious to do, but why isn't there an observafion-based analysis in the manuscript to\ncomplement the modeling results (e.g. using the Dai data set)? Related to this issue, the text should be\nadapted to insert, in all instances, the word \"modeled\" when referring to snowfall probabilifies and\nother model-derived quanfifies.\n2) Throughout the paper there seems to be an over-use of abbreviafions which make the paper\nunnecessarily challenging to read - the secfion starfing on line 190 is a prime example.\n3) Greater emphasis needs to be placed on the analysis of Snow Dominated Events for situafions with\npre-exisfing snow cover versus events without pre-exisfing snow cover. Wouldn't the presence of snow\nlargely negate land cover associated differences between urban versus buffer areas as everything is\ncovered in highly reflecfive snow?\n4) What is the impact of the background land cover type of the buffer areas? Given that albedo\ndifferences between pair members seems to play a large role in determining differences in snowfall\nprobability, it would stand to reason that the albedo of buffer areas would be important. For example, a\nforest typically has a lower albedo than a grassland. In addifion, a small amount of snow cover will\ndramafically alter the albedo grassland whereas it will have a more subtle impact on the albedo of a\nforest. These differences in buffer albedo would then have an appreciable impact on the difference\nbetween the buffer albedo and an adjacent urban area. Can this be evaluated with the data presented\nherein? Can the analysis be strafified to look at the differences in snowfall probability across differences\nin background buffer land cover type? Or at the very least, this effect seems worthy of discussion within\nthe manuscript as the spafial pafterns in daSP seem to evoke a land cover control.\n5) In the secfion starfing on line 190 related to trends: Can you report the trend in daSP? Is it possible\nthat if warming potenfially pushes more pairs toward the rain-snow threshold temperature, might the\nsubtle differences in temperature between pair members cause a greater decrease in snowfall\nprobability an urban areas relafive to the corresponding buffer area? I think the analysis suggests there is\ngreater sensifivity in the urban areas but an explicit analysis of daSP trend seems warranted."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "The manuscript has been improved significantly in the current version. Specially, the authors added\ncomplementary materials and explanafions to the data and methods part. Therefore, the manuscript is\nclear and readable. The findings are worth publishing in my understanding based on its novelty and solid\nanalysis. However, I sfill have some suggesfions.\n(1) The keywords should befter contain the informafion of “urbanizafion”.\n(2) Numerical values should be added in the abstract to make the results more interesfing.\n(3) The conclusion is not impressive enough due to its length.\n(4) Detail of text often requires aftenfion.", "reviewer_3": "MAJOR COMMENT 1: Importance of the analysis.\nWhy is this work important? Why is it important to the public to know that urbanizafion impacts\nsnowfall fracfion? Reviewer 1 raised this quesfion and I agree that the authors need to make the case as\nto why this work is important. Urban areas are not water sources so the snowfall/rainfall disfincfion is\narguably not important in a water resources context. In my view, it is perhaps most important in terms of\nhow the public experiences the change in precipitafion type and the degree to which it is influenced by\nland cover change versus climate change.\nFor example, the final sentence of the abstract states:\n“Results highlight a clear imprint of urbanizafion on precipitafion phase and underscore the need to\nconsider these influences while projecfing hydro-meteorological risks.”\nWhat do the authors mean by “projecfing hydro-meteorological risks”? I don’t know what that means\nand it seems like an overtly vague statement. While this one sentence is a very specific issue, it is\nperhaps the most important sentence in the paper and begs the quesfion as to why this work is\nimportant beyond random curiosity.\nMAJOR COMMENT 2: Misleading language about what is “observed” versus “modeled”.\nAs with my comment on the previous review, the paper needs to be explicit that the results are based on\na model and not observafions. The authors responded that they have clarified this but I do not agree\nthat it has been adequately addressed as nowhere does the word “model” or “modeled” appear within,\nor prior to, the text reporfing the results. And in fact, the authors use the word(s) “observed” or “were\nobserved” when reporfing the results but in fact, these are “modeled” not “observed” quanfifies.\nMAJOR COMMENT 3: Improving the robustness of the analyses by ufilizing observed data.\nThe point I raised in the inifial review about the lack of observafion-based analyses has not been\nadequately addressed by the authors.\nMy inifial comment was: “The authors rely on a model and not observafions. I assume there is a good\njusfificafion for this as it seems obvious to do, but why isn't there an observafion-based analysis in the\nmanuscript to complement the modeling results (e.g. using the Dai data set)? Related to this issue, the\ntext should be adapted to insert, in all instances, the word \"modeled\" when referring to snowfall\nprobabilifies and other model-derived quanfifies.”\nThe authors responded as follows:\n“Of the total 7,415 UBPs, while 1,973 observafion points fall within the urban regions and 246 in buffer\nregions, unfortunately there are zero UBPs with at least one stafion within both urban and\ncorresponding buffer region. Therefore, the dataset cannot be used for the analysis carried out in this\nstudy.”"}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have appropriately addressed my final comment and I believe the paper can be published\nwithout further revision - however, I strongly urge the authors and the editorial staff to double check the\npaper for typos and / or grammafical errors (as a standard pracfice)."}, "author_response": ""}]

Decision letter