The Marcellus Shale is considered to be the largest unconventional shale-gas resource in the United States. Two critical factors for unconventional shale reservoirs are the response of a unit to hydraulic fracture stimulation and gas content. The fracture attributes reflect the geomechanical properties of the rocks, which are partly related to rock mineralogy. The natural gas content of a shale reservoir rock is strongly linked to organic matter content, measured by total organic carbon (TOC). A mudstone lithofacies is a vertically and laterally continuous zone with similar mineral composition, rock geomechanical properties, and TOC content. Core, log, and seismic data were used to build a three-dimensional (3-D) mudrock lithofacies model from core to wells and, finally, to regional scale. An artificial neural network was used for lithofacies prediction. Eight petrophysical parameters derived from conventional logs were determined as critical inputs. Advanced logs, such as pulsed neutron spectroscopy, with log-determined mineral composition and TOC data were used to improve and confirm the quantitative relationship between conventional logs and lithofacies. Sequential indicator simulation performed well for 3-D modeling of Marcellus Shale lithofacies. The interplay of dilution by terrigenous detritus, organic matter productivity, and organic matter preservation and decomposition affected the distribution of Marcellus Shale lithofacies distribution, which may be attributed to water depth and the distance to shoreline. The trend of normalized average gas production rate from horizontal wells supported our approach to modeling Marcellus Shale lithofacies. The proposed 3-D modeling approach may be helpful for optimizing the design of horizontal well trajectories and hydraulic fracture stimulation strategies.