This study presents a hybrid framework to predict stability solutions of buried structures under active trapdoor conditions in natural clays with anisotropy and heterogeneity by combining physics-based and data-driven modeling. Finite-element limit analysis (FELA) with a newly developed anisotropic undrained shear (AUS) failure criterion is used to identify the underlying active failure mechanisms as well as to develop a numerical (physics-based) database of stability numbers for both planar and circular trapdoors. Practical considerations are given for natural clays to three linearly increasing shear strengths in compression, extension, and direct simple shear in the AUS material model. The obtained numerical solutions are compared and validated with published solutions in the literature. A multivariate adaptive regression splines (MARS) algorithm is further utilized to learn the numerical solutions to act as fast FELA data-driven surrogates for stability evaluation. The current MARS-based modeling provides both relative importance index and accurate design equations that can be used with confidence by practitioners.