The greatest challenges of rigorously modeling coupled hydro-mechanical (HM) processes in fractured geological media at different scales are associated with computational geometry. These challenges include dynamic shearing and opening of intersecting fractures at discrete fracture scales as a result of coupled processes, and contact alteration along rough fracture surfaces that triggers structural and physical changes of fractures at micro-asperity scale. In this paper, these challenges are tackled by developing a comprehensive modeling approach for coupled processes in fractured geological media based on numerical manifold method (NMM) at multiple scales. Based on their distinct geometric features, fractures are categorized into three different scales: dominant fracture, discrete fracture, and discontinuum asperity scales. Here the scale is relative, that of the fracture relative to that of the research interest or domain. Different geometric representations of fractures at different scales are used, and different governing equations and constitutive relationships are applied. For dominant fractures, a finite thickness zone model is developed to treat a fracture as a porous nonlinear domain. Nonlinear fracture mechanical behavior is accurately modeled with an implicit approach based on strain energy. For discrete fractures, a zero-dimensional model was developed for analyzing fluid flow and mechanics in fractures that are geometrically treated as boundaries of the rock matrix. With the zero-dimensional model, these fractures can be modeled with arbitrary orientations and intersections. They can be fluid conduits or seals, and can be open, bonded or sliding. For the discontinuum asperity scale, the geometry of rough fracture surfaces is explicitly represented and contacts involving dynamic alteration of contacts among asperities are rigorously calculated. Using this approach, fracture alteration caused by deformation, re-arrangement and sliding of rough surfaces can be captured. Our comprehensive model is able to handle the computational challenges with accurate representation of intersections and shearing of fractures at the discrete fracture scale and rigorously treats contacts along rough fracture surfaces at the discontinuum asperity scale. With future development of three-dimensional (3D) geometric representation of discrete fracture networks in porous rock and contacts among multi-body systems, this model is promising as a basis of 3D fully coupled analysis of fractures at multiple scales, for advancing understanding and optimizing energy recovery and storage in fractured geological media.