A fourth-order in space and second-order in time 3D staggered (SG) and rotated-staggered-grid (RSG) method for the solution of Biot's equation are presented. The numerical dispersion and stability conditions are derived using a von Neumann analysis. The exact stability condition is calculated from the roots of a 12th-order polynomial and therefore no nontrivial expression exists. To overcome this, a 1D stability condition is usually generalized to three dimensions. It is shown that in certain cases, the 1D approximate stability condition is violated by a 3D SG method. The RSG method obeys the approximate 1D stability condition for the material properties and spatiotemporal scales in the examples shown. Both methods have been verified against an analytical solution for an infinite homogeneous porous medium with a misfit error of less than 0.5%. A free surface has been implement-ed to test the accuracy of this boundary condition. It also serves as a test of the methods to include high material contrasts. The methods have been compared with a quasi-analytical solution. For the specific material properties, spatial grid scaling, and propagation distance used in the test, a maximum error of 3.5% for the SG and 4.1% for the RSG was found. These errors depend on the propagation distance, temporal and spatial scales, and accuracy of the quasi-analytical solution. No discernable difference was found between the two methods except for time steps comparable with the stability-criteria threshold time step, the SG was found to be unstable. However, the RSG remained stable for a homogeneous half-space. Time steps, comparable to the stability criteria, reduce the computational time at the cost of a reduction in accuracy. The methods allow wave propagation to be modeled in a porous medium with a free surface.