Many geological materials are multi-scale in nature. They typically demonstrate a nonlinear and path-dependent behavior, which is determined by the microscopic heterogeneities as well as properties, shape, and distribution of the constituents. One of the major challenges associated with multi-scale modeling of complex geological systems includes computational costs of high-fidelity simulations and transferring the information across different length scales. Consequently, design and analysis of these systems often involves phenomenological continuum-scales constitutive models that are implemented in a computational framework such as the Finite Element Method (FEM). Recently, data- driven techniques have emerged as a promising alternative to the traditional constitutive models, which can expedite simulations of complex geological systems. In this work, we present the recent advances in machine learning-enabled surrogate modeling for nonlinear and path-dependent behavior of geomaterials and their application in the computational framework of FEM. Simulation results show that the proposed model is able to capture the path-dependent behavior of brittle and/or ductile geomaterials. Numerical results are presented to study the effects of model architecture, training parameters and training datasets on the accuracy of FEM simulations.