Abstract The Poisson-Boltzmann (PB) model is an effective approach for the electrostatics analysis of solvated biomolecules. The nonlinearity associated with the PB equation is critical when the underlying electrostatic potential is strong, but is extremely difficult to solve numerically. In this paper, we construct two operator splitting alternating direction implicit (ADI) schemes to efficiently and stably solve the nonlinear PB equation in a pseudo-transient continuation approach. The operator splitting framework enables an analytical integration of the nonlinear term that suppresses the nonlinear instability. A standard finite difference scheme weighted by piecewise dielectric constants varying across the molecular surface is employed to discretize the nonhomogeneous diffusion term of the nonlinear PB equation, and yields tridiagonal matrices in the Douglas and Douglas-Rachford type ADI schemes. The proposed time splitting ADI schemes are different from all existing pseudo-transient continuation approaches for solving the classical nonlinear PB equation in the sense that they are fully implicit. In a numerical benchmark example, the steady state solutions of the fully-implicit ADI schemes based on different initial values all converge to the time invariant analytical solution, while those of the explicit Euler and semi-implicit ADI schemes blow up when the magnitude of the initial solution is large. For the solvation analysis in applications to real biomolecules with various sizes, the time stability of the proposed ADI schemes can be maintained even using very large time increments, demonstrating the efficiency and stability of the present methods for biomolecular simulation.