ABSTRACT: Injection-induced seismicity (IIS) depends on pore pressure, in-situ stress state, and fault orientation; generally occurs in basement rock that contains fractures and faults; and moves away from the injection well as a nonlinear diffusion process. Therefore, to numerically model IIS a code should incorporate flow and geomechanics, the presence of fractures and faults, and the capability for hydraulic diffusivity to evolve with effective stress and failure history. In this work, we introduce and verify a modeling framework that allows hydraulic diffusivity to evolve as fractures open and close. Details and challenges in code development are discussed, including how the Bandis model for normal fracture deformation can be used to calculate hydraulic diffusivity as a function of effective normal stress. The discrete fracture network and matrix (DFNM) model is implemented in PFLOTRAN such that hydraulic diffusivity has different constitutive relationships for fracture and matrix grid cells. This model is applied to understand the recent IIS near Greeley, Colorado, and its results are compared to: (a) a traditional DFNM model where hydraulic diffusivity cannot evolve and (b) an equivalent porous media (EPM) model where the effect of the fractures are averaged over a large region of rock. The new DFNM model predicts critical pressure will propagate farther from an injection well. This modeling framework shows promise for applications where fracture and matrix flow are important and hydraulic diffusivity is a function of pressure, stress, and/or shear failure history.