This study presents comprehensive coupled hydro-mechanical numerical modeling of the first and second hydraulic stimulations at the Pohang enhanced geothermal system (EGS) site in order to improve the understanding on the key stimulation mechanisms in the fractured reservoir. Two models for PX-2 and PX-1 wells were developed including fracture zones of different orientations and permeability. Shear dilation and frictional plastic strain-softening were implemented in a non-linear stress dependent fracture model. A history matching of the wellhead pressure curves was carried out in the early days of each stimulation. The PX-2 model with a fracture jacking mechanism achieved the reproduction of highly reversible non-linear aperture changes observed in situ. The PX-1 model successfully simulated a wellhead pressure peak and drop associated with shear slip and dilation by applying frictional plastic strain softening. Simulated critical wellhead pressure for shear slip near the PX-1 well was greater than the prediction by a simple slip potential analysis due to local stress changes by the poroelastic effect. The numerical modeling confirmed that the combination of shear dilation and jacking was adequate to capture the pressure evolution during the increasing step-rate test at PX-1. Besides, spatio-temporal changes in pressure, total stress, and permeability evaluated in the fracture zones greatly enhanced the understanding on the coupled behavior caused by the hydraulic stimulations. The possibility of zones with lower permeability away from the PX-1 well was suggested by an alternative model. The current numerical study demonstrates that the key hydro-mechanical processes of shear slip and dilation and hydraulic jacking observed in the fractured reservoir can be successfully reproduced.