<p>Accurately simulating the complex response of foundations in geotechnical engineering relies heavily on the numerical implementation of elastoplastic constitutive models for natural clays. This paper presents a numerical realization format for the elastic–plastic constitutive model of the yield surface under both upper and lower loading conditions, incorporating the plastic anisotropy characteristics, structural properties, and evolution laws inherent to natural soft clay. An automatic error control sub-incremental step display integration algorithm is introduced. A user defined material (UMAT) subroutine was developed utilizing the sub-incremental step display integration algorithm for automatic error control, and it has been successfully integrated into the finite element software ABAQUS. Through the triaxial simulation test, the results obtained from the model were compared with the triaxial test data. The findings confirm a high degree of agreement between the predicted outcomes and experimental measurements, thereby validating the robustness of the model. Furthermore, when applied to a practical case involving embankment filling, the predictions generated by this model exhibit strong consistency with measured data. This further substantiates its practical application value and highlights its potential in addressing complex boundary value problems within geotechnical engineering. It is noteworthy that factors such as soil initial anisotropy, structural characteristics, overconsolidation, and their evolving parameters exert varying degrees of influence on surface settlement and soil pore water pressure.</p>