In the present study a numerical model developed by Lynett and Liu (2002) is modified to include density difference in a stratified two-layer fluid in a three-dimensional internal wave domain. The internal solitary wave (ISW) in the model is assumed to be weakly nonlinear and weakly dispersive, and the viscosity effects at all boundaries are ignored. The governing equations based on the Navier-Stokes and Euler equations are solved for internal solitary wave propagation over variable seabed topography. Theoretical formulations are established, from which analytical solutions are obtained, in addition to numerical results. Wave profiles from previous experimental studies are compared with the numerical results from the present analytical solutions. Numerical models developed on the basis of the present analytical solutions are better than those developed by Lynett and Liu (2002). The results of numerical modeling agree well with the experimental data. In the present study a numerical model developed by Lynett and Liu (2002) is modified to include density difference in a stratified two-layer fluid in a three-dimensional internal wave domain. The internal solitary wave (ISW) in the model is assumed to be weakly nonlinear and weakly dispersive, and the viscosity effects at all boundaries are ignored. The governing equations based on the Navier-Stokes and Euler equations are solved for internal solitary wave propagation over variable seabed topography. Theoretical formulations are from, from which analytical solutions are produced, in addition to numerical results from the present analytical solutions. The numerical models developed on the basis of the present analytical solutions are better than those developed by Lynett and Liu (2002) The results of numerical modeling agree well with the experimental data.