This work presents a methodology for estimating seasonal snow water equivalent (SWE) from the use of remotely sensed Visible and Near Infrared observations from the Landsat mission. The method is comprised of two main components: (1) a coupled land surface model and snow depletion curve model, which is used to generate an ensemble of predictions of SWE and snow cover area for a given set of (uncertain) inputs, and (2) a reanalysis step, which updates estimation variables to be consistent with the satellite observed depletion of the fractional snow cover time series. This method was applied over the Sierra Nevada (USA) based on the assimilation of remotely sensed fractional snow covered area data over the Landsat 5-8 record (1985-2016). The verified dataset (based on a comparison with over 9000 station years of in situ data) exhibited mean and root-mean-square errors less than 3 and 13 cm, respectively, and correlations with in situ SWE observations of greater than 0.95. The method (fully Bayesian), resolution (daily, 90-meter), temporal extent (32 years), and accuracy provide a unique dataset for investigating snow processes. In particular, this presentation illustrates how the reanalysis dataset was used to provide climatology of the seasonal snowfall accumulation rates, distributions, and variability over the last three decades.