We present a novel three dielectric layer hybrid solvation model for treating electrostatic interactions of biomolecules in solvents using the Poisson-Boltzmann equation. In this model, the interior spherical cavity contains the solute and explicit solvent molecules. An intermediate buffer layer is introduced, which also contains solvent molecules. Outside the spherical shell defines the exterior layer, where bulk solvent is modeled implicitly and characterized by a dielectric constant. Within the intermediate layer, a special dielectric permittivity profile is constructed to give a continuous transition from the interior cavity to the exterior layer. The selection of this special profile using a harmonic interpolation allows an analytical solution of the model by generalizing the classical Kirkwood series expansion. To speed up numerical calculations of the electrostatic potential solutions, discrete image charges are employed following previous work [1]. Two approaches for constructing discrete image approximations to the potentials are considered: Semi-analytical and least square methods. Both methods are employed for the reaction field of solvents without and with finite ionic strength. Numerical results are presented to validate the accuracy and effectiveness of the image charge methods. This work is supported by NIH 1R01 GM083600-02. Z. Xu is also partially supported by the Charlotte Research Institute through a Duke Postdoctoral Fellowship.[1] W. Cai, S. Deng and D. Jacobs, J. Comp. Phys. 223, 846-864 (2007).