High-resolution earthquake simulations are essential for understanding seismic wave propagation in complex geological structures. However, finite element simulations require significant computational resources, particularly when exploring multiple scenarios or accounting for uncertainties. This PhD project aims to develop causal surrogate models that reduce computational costs while preserving the physical consistency of seismic simulations. The proposed approach combines causal machine learning, Bayesian statistics, and high-performance computing (HPC). Using data generated by large-scale finite element simulations, the research will investigate methods for identifying directed relationships between seismic variables across space and time. These causal relationships will then be used to construct surrogate models capable of approximating selected components of numerical simulators. Bayesian approaches will also be explored to quantify uncertainties in the learned causal structures and model predictions. The developed methods will be integrated into parallel simulation environments, particularly ArcaneFEM and PSD, to ensure their applicability to large-scale problems. The expected results include efficient causal learning algorithms, computationally scalable surrogate models, and hybrid simulation strategies combining numerical methods and artificial intelligence. Ultimately, this research aims to improve the efficiency of earthquake simulations and support uncertainty quantification, adaptive mesh refinement, and data assimilation.