This work focuses on the numerical modeling and optimization of thermal energy storage. The system proposed in this study consists of phase change materials (PCMs) with spherical heat transfer enhancers. In order to optimize the position and size of the enhancers, multiple numerical simulations must be performed. However, it is well known that numerical simulations of phase change heat transfer processes in porous media are computationally expensive and time-consuming [1, 2]. Therefore, in order to reduce the computational cost associated with solving the heat conduction equations that take phase change into account, we propose a efficient reduced-order model. More specifically, the reduced basis is generated using the Proper Orthogonal Decomposition (POD) applied to a set of snapshots obtained from the full-order model at selected time instants. To handle the nonlinear terms in the governing equations, several strategies are investigated, including standard POD-based approaches, the Discrete Empirical Interpolation method (DEIM), and neural network-based techniques. Finally, the accuracy and computational efficiency of the proposed reduced-order models are validated through multiple numerical examples. Future work will focus on extending the proposed framework to reduced-order modeling of parameterized versions.