The hydrological effects of earthquakes are determined by the style of fault displacement, rather than simply by the magnitude of the earthquake. In the past, many researchers used the analytical solution of Okada (1992) to estimate the pore‐pressure field, which is derived through the stress field by Skempton’s B coefficient, after the stress field is calculated from a given forced fault slip in Okada’s method (Okada, 1992). This approach is a one‐way coupling approach, as fluid has no effect on rock behavior, and Skempton’s B coefficient must be given separately to know the pore‐pressure field. We create a fault‐slip model in conjunction with a poroelastic model and present a new approach in which fault movement is a consequence of reduced friction occurring in a saturated porous medium. This model represents a numerical parallel to the traditional Okada’s approach, but with two advantages: (1) it is well suited for more complex geometries and heterogeneous formations for which analytical methods fail; (2) it is fully coupled so that fluid effects are also taken into account. This model shows that all three different faulting regimes exhibit completely different pore‐pressure change distributions in 3D space. The transient pore‐pressure change after a strike‐slip faulting event is axis symmetrical with respect to the imaginary vertical line passing through the epicenter and hypocenter, with increased and decreased areas evenly distributed in a juxtaposed pattern. The transient pore‐pressure changes associated with normal and thrust faulting all display a circular pattern on the model surface, with the largest pore‐pressure increase occurring at the epicenter for normal faulting and the largest decrease for thrust faulting. Many observed field phenomena have been explained in surprising detail by this model.