Liquid metal blankets play a critical role in fusion reactor design, fulfilling functions such as tritium breeding, plasma energy conversion into heat, and shielding sensitive reactor components from harmful radiation. Designing and optimizing these blankets necessitate extensive magnetohydrodynamic (MHD) simulations that capture the complex interaction between liquid metal flows and the strong magnetic fields used to confine the plasma. However, the computational cost of such simulations is substantial even on advanced supercomputers, making it challenging to address optimization and uncertainty quantification tasks that require a large number of forward model evaluations. Surrogate modeling provides a practical solution by approximating the expensive model with computationally cheaper alternatives constructed from a limited set of high-fidelity simulations. We investigate several popular surrogate modeling approaches for MHD flows in liquid metal fusion blankets, including sparse grid, sparse polynomial expansion, and Gaussian process regression. Focusing on MHD flow simulations within a fixed square domain, we consider Hartmann numbers ranging from 100 to 10,000 and wall conductivities from 1.E-5 to 1. Results demonstrate that all three methods can achieve relative errors below 0.1%, using 200-600 training simulations. In our study, sparse grid interpolation delivers the highest accuracy per simulation, while sparse polynomial expansion and Gaussian process regression offer greater flexibility in sampling input selection at the expense of minor accuracy loss. Additionally, preliminary tests with neural network-based approaches reveal that these methods require substantially larger datasets and provide no tangible benefits when simulations are limited in number.