In this paper, three methods are used in order to obtain the solution for the propagation of water solitons over finite and variable depth. First, the exact analytical solitary wave solutions of the one-dimensional non-linear Boussinesq equations under shallow water condition are described for constant and variable depth. Second, the three-dimensional Fully Non-linear Potential Flow code OceanWave3D is used in order to obtain the numerical solutions for the solitary waves’ propagation over same depth ranges, providing robust solutions for the potential flow problem. Third, Computational Fluid Dynamics’ tool OpenFOAM is used in order to obtain the viscous solution for the same problem, however, without the accounts of turbulence models. The free-surface profiles are drawn and compared; and the stability of the numerical solutions are assessed. Since the approximations of Boussinesq-type equations depend mainly on the orders of magnitude of amplitude and depth, the numerical-analytical comparison will draw the limits for the validity of the analytical solutions. On the other hand, the comparison will provide the limits where viscous effects start playing an important role, whereas the CFD simulations predict the occurrence of wave breaking. These benchmark cases are compared with past references. After all, results regarding the same phenomena have been described in the literature according to, e.g. Fully Non-linear Boussinesq Models, and Fully Nonlinear Potential Flow schemes solved by Boundary Element Methods. Last but not least, the open source Fully Non-linear Potential Flow code is used in order to provide the potential flow solution for some extra cases of water soliton propagation, in order to capture the trends in weak shoaling scenarios.