Snow on Arctic sea ice plays many roles in Arctic climate feedbacks; in particular, through its influence on the underlying sea ice. Estimates of snow depth on Arctic sea ice are also a key input for deriving sea ice thickness from satellite lidar altimetry measurements, such as those from ICESat-2. Snow depth throughout the Arctic basin is difficult to observationally quantify due to the limited extent of in-situ measurements and biases in remote observations. Snow-on-sea-ice models, such as the NASA Eulerian Snow On Sea Ice Model (NESOSIM) can produce basin-wide estimates of snow depth and density on Arctic sea ice. NESOSIM version 1.1 is a 2-layer model with simple representations of snow accumulation, wind packing, loss due to blowing snow, and redistribution due to sea ice motion. Reanalysis snowfall input to NESOSIM is scaled to observed snowfall derived from CloudSat satellite radar measurements. We use a Metropolis Markov Chain Monte Carlo approach to indirectly calibrate NESOSIM model free parameters to snow depth observations from Operation IceBridge and CRREL-Dartmouth snow buoys, and historical density measurements from Soviet drifting stations. This approach produces estimates of the parameters and their probability distributions, from which the contribution of model parameter uncertainty to snow depth uncertainty can be estimated. We further incorporate an estimate of uncertainty due to model sensitivity to the choice of snowfall input, and construct a blended snow-on-sea-ice product that combines the output from NESOSIM driven with ERA5, MERRA-2, and JRA55 snowfall products. Finally, we briefly examine the impact of this calibration on sea ice thickness derived using NESOSIM output and ICESat-2 freeboard measurements.