||The energy industry is facing major challenges as world energy demand continues to increase. Unconventional resources, including oil shale, have the capacity to meet that demand as production technology develops. There are many challenges associated with production of fuels from these resources including economics, environmental considerations, sustainability, and government policies. Improving predictive capabilities for energy production from these resources could play a major role in addressing these challenges. Simulation of thermal/reactive reservoir systems is complex. Heat transfers through the rock and fluids by conduction and convection, chemical reactions occur, phases change, multiphase fluids ow, and rock mechanical properties change. Each of these physical processes occurs simultaneously, affecting each other. Important physical processes also occur at widely varying length and time scales. Even where physical process models may be appropriate, the input data are highly uncertain. Useful modeling tools must find a balance between solution accuracy and computational effeciency. Simplifying assumptions are typically made according to data and experience. However, these simplifying assumptions may not be justified where data are sparse and where system characteristics change as a process unfolds. This research explores methods to expose the most important physical parameters and models for making efficient and useful predictions. Reservoir simulation methods for approaching thermal and reactive problems have been explored and analyzed showing heating by conduction alone is slow, and removal of products in a rubble system is not trivial. Experimental design methods have been implemented to expose some important parameters for predicting liquid fuel production, and surrogate models have been created with those parameters that approximate full simulation within 15% accuracy. Expected variations in kinetics and relative permeability modeling parameters predict a normal probability with a mean of 295 bbls (about 39 wt% initial kerogen in place) ultimate oil recovery with a standard deviation of 40 bbls (about 5 wt% initial kerogen in place). A case study involving a hybrid process involving in situ pyrolysis, in situ combustion, and CO2 enhanced oil recovery has shown that the energy input requirement for oil shale heating can be reduced by more than 335 million BTU, or 55%, and produce 160 bbls or 36% more oil. Finally, novel methods for performing these complex simulations are discussed. Progress and challenges with these novel methods are also discussed.