The modeling of stress-strain behavior of geomaterials, such as soils, is key to the accurate analyses of complicated geotechnical engineering structures. Traditional elastoplastic modeling concepts, characterized by a single yield surface, however, limit our ability to model complex stress-strain responses. In this dissertation, a novel modeling concept called the Middle Surface Concept (MSC), is developed using multiple pseudo yield surfaces. The MSC is first developed to model saturated sand behavior under monotonic triaxial loading conditions and then extended to the general stress space. Single element model predictions are compared to laboratory tests results for three different types of sands subjected to various loading conditions and reasonable comparisons are obtained. In order to implement the general stress space MSC sand model into a finite element method, the consistent tangent stiffness matrix is developed and the model is numerically integrated using the generalized trapezoidal rule. Some useful restrictions in terms of Poisson's ratio for various flow rules used in constitutive models for granular materials are also developed. The MSC sand model is implemented into a fully coupled computer code, DYSAC2, and predictions are made for a centrifuge model subjected to base shaking. Reasonable comparisons between DYSAC2 predictions and centrifuge model test results are obtained validating the performance of the MSC sand model in boundary value problems. Finally, the MSC is expanded to model unsaturated sand or silt behavior under triaxial monotonic loading conditions. Two pseudo yield surfaces are utilized to model the effects of suction on the stress-strain behavior of unsaturated sands and silts.