  In our previous work [SIAM J. Sci. Comput. 43(3) (2021) B784-B810], an
accurate hyper-singular boundary integral equation method for dynamic
poroelasticity in two dimensions has been developed. This work is devoted to
studying the more complex and difficult three-dimensional problems with Neumann
boundary condition and both the direct and indirect methods are adopted to
construct combined boundary integral equations. The strongly-singular and
hyper-singular integral operators are reformulated into compositions of
weakly-singular integral operators and tangential-derivative operators, which
allow us to prove the jump relations associated with the poroelastic layer
potentials and boundary integral operators in a simple manner. Relying on both
the investigated spectral properties of the strongly-singular operators, which
indicate that the corresponding eigenvalues accumulate at three points whose
values are only dependent on two Lam\'e constants, and the spectral properties
of the Calder\'on relations of the poroelasticity, we propose
low-GMRES-iteration regularized integral equations. Numerical examples are
presented to demonstrate the accuracy and efficiency of the proposed
methodology by means of a Chebyshev-based rectangular-polar solver.
