The spectral distribution of marine remote sensing reflectance, Rrs, is the fundamental measurement of ocean color science, from which a host of bio-optical and biogeochemical properties of the water column can be derived. Estimation of uncertainty in these derived properties is thus dependent on knowledge of the uncertainty in satellite-retrieved Rrs (uc(Rrs)) at each pixel. Uncertainty in Rrs, in turn, is dependent on the propagation of various uncertainty sources through the Rrs retrieval process, namely the atmospheric correction (AC). A derivative-based method for uncertainty propagation is established here to calculate the pixel-level uncertainty in Rrs, as retrieved using NASA's multiple-scattering epsilon (MSEPS) AC algorithm and verified using Monte Carlo (MC) analysis. The approach is then applied to measurements from the Moderate Resolution Imaging Spectroradiometer (MODIS) on the Aqua satellite, with uncertainty sources including instrument random noise, instrument systematic uncertainty, and forward model uncertainty. The uc(Rrs) is verified by comparison with statistical analysis of coincident retrievals from MODIS and in situ Rrs measurements, and our approach performs well in most cases. Based on analysis of an example 8-day global products, we also show that relative uncertainty in Rrs at blue bands has a similar spatial pattern to the derived concentration of the phytoplankton pigment chlorophyll-a (chl-a), and around 7.3%, 17.0%, and 35.2% of all clear water pixels (chl-a ≤ 0.1 mg/m3) with valid uc(Rrs) have a relative uncertainty ≤ 5% at bands 412 nm, 443 nm, and 488 nm respectively, which is a common goal of ocean color retrievals for clear waters. While the analysis shows that uc(Rrs) calculated from our derivative-based method is reasonable, some issues need further investigation, including improved knowledge of forward model uncertainty and systematic uncertainty in instrument calibration.