The safe operation of nuclear reactors during severe accidents relies on the effective cooling of debris beds. Debris beds consist of fragmented and molten core materials that accumulate following significant core damage. Existing safety assessment models typically assume uniform debris bed structures. However, experimental and observational evidence indicates that debris beds often exhibit internal heterogeneities, such as regions with reduced porosity known as cake layers. These heterogeneities can significantly influence coolant flow patterns and quenching behavior. This thesis addresses the challenge of incorporating such structural variations into predictive models. A hybrid framework is developed, combining high-fidelity mechanistic simulations of two-phase flow and heat transfer in debris beds with machine learning-based surrogate modeling. The simulations explicitly account for spatially varying porosity and the presence of cake layers. Data generated from these simulations are used to train surrogate models capable of rapidly predicting quenching performance under diverse structural conditions. This approach enables substantial reductions in computational expense while maintaining physical accuracy. The objective is to establish a scalable and physically consistent methodology for evaluating debris bed coolability in the presence of internal heterogeneities. This work seeks to enhance understanding of the effects of porosity stratification on coolant transport and dryout risk. Additionally, it aims to develop efficient surrogate tools suitable for parametric studies and real-time safety evaluations. By integrating detailed physics-based simulations with data-driven modeling, this research contributes to advancing nuclear accident management strategies. Furthermore, it supports the deployment of digital twin technologies for improved reactor safety monitoring and design optimization.