Geostress, the stress within the Earth’s crust, is a critical factor in subsurface systems, yet its direct measurement is challenging and often prohibitively costly. This study develops a mathematical framework to model geostress distribution trends under quasi-static conditions using finite element methods (FEM). Leveraging known surface deformation, subsurface horizon deformation, and fault geometry—readily obtainable via InSAR, 3D seismic interpretation, or other geophysical techniques—we derive a model based on gravitational loading and gravity-derived prestress, yielding a standardized trend for underground horizon geostress distribution. Mathematical derivations address regional variations in material properties, anisotropy, layer attributes, and fault parameters, enhancing the framework’s theoretical flexibility. Dynamic factors, such as tectonic events and injection-induced seismicity, challenge this static model, exposing its limitations and raising significant mathematical questions for future theoretical exploration. These findings establish a robust foundation for advancing geostress analysis, offering insights into subsurface stress modeling through rigorous mathematical methods.