Thermal energy storage (TES) modules are specifically designed to respond to transient thermal loading. Their dynamic response depends on the overall structure of the module, including module geometry and dimensions, the internal spatial distribution of phase change material (PCM) and conductive heat-spreading elements, and the thermophysical properties of the different materials composing the module. However, due to the complexity of analyzing a system’s dynamic thermal response to transient input signals, optimal design of a TES module for a particular application is challenging. Conventional design approaches are limited by (1) the computational cost associated with high fidelity simulation of heat transfer in nonlinear systems undergoing a phase transition and (2) the lack of model integration with robust optimization tools. To overcome these challenges, I derive reduced-order dynamic models of two different metal-PCM composite TES modules and validate them against a high fidelity CFD model. Through simulation and validation of both turbulent and laminar flow cases, I demonstrate the accuracy of the reduced-order models in predicting, both spatially and temporally, the evolution of the dynamic model states and other system variables of interest, such as PCM melt fraction. The validated models are used to conduct univariate and bivariate parametric studies to understand the effects of various design parameters on different performance metrics. Finally, a case study is presented in which the models are used to conduct detailed design optimization for the two HX geometries.