A Bayesian approach to estimate bias and uncertainty in snowfall precipitation from MERRA-2 and other precipitation products was applied over High Mountain Asia (HMA), using a newly developed snow reanalysis method. Starting from an ‘uninformed’ prior probability distribution, a posterior scaling factor applied to MERRA-2 snowfall was derived by constraining model-based estimates of seasonal snow accumulation and ablation over the water year (WY) with fractional snow covered area (fSCA) measurements derived from Landsat and MODIS (MODSCAG). Several sub-domains (nine representative 1° by 1° tiles) across HMA were examined over the period WYs 2001-2015 and compiled into an uncertainty parameterization where a lognormal distribution was fitted to the empirical posterior distribution with a mean of 1.54 (median of 1.19) and coefficient of variation of 0.83, indicating that MERRA-2 underestimates snowfall on average by ~54% with sizeable uncertainty. For reference, the uncertainties in snowfall precipitation from the ERA5 and APHRODITE-2 precipitation produces were also evaluated, and these products were found to underestimate snowfall, on average by a factor around 1.78 and 3.34 (with median scaling factors of 1.42 and 2.51) respectively. The results indicate that snowfall precipitation at high-elevations dominated by snowfall is underestimated in most existing products, especially in the gauge-based APHRODITE-2 product, where the biases were also found to exhibit geographical variations with the largest underestimation in monsoon-influenced high-elevation tiles. The derived MERRA-2 uncertainty model is being used to develop a full domain-wide HMA snow reanalysis, which will shed further light onto the space-time variations in snowfall biases in these products.