This paper presents the development of a theoretical model of fully nonlinear and weakly dispersive (FNWD) waves and numerical techniques for simulating the propagation, interaction, and transformation of solitary waves. Using the standard expansion method and without the limit of small nonlinear parameter defined as the ratio of the wave height versus water depth, a set of model equations describing the FNWD waves in a domain of moderately varying bottom topography are formulated. Exact solitary wave solutions satisfying the FNWD equations are also derived. Numerically, a time-accurate and stabilized finite-element code to solve the governing equations is developed for wave simulations. The solitary wave solutions of FNWD, weakly nonlinear and weakly dispersive (WNWD), and Laplace equations based models in terms of wave profile and phase speed are compared to examine their related features and differences. Investigations on the overtaking collision of two unidirectional solitary waves of different amplitudes, i.e.,1α and2α where12α>α, are carried out using both the FNWD and WNWD water wave models. Selected cases by running the FNWD and WNWD models are performed to identify the critical values of12α/α for forming a flattened merging wave peak, which is the condition used to determine if the stronger wave is to pass through the weaker one or both waves are to remain separated during the encountering process. It is interesting to note the critical values of12α/α obtained from the FNWD and WNWD models are found to be different and greater than the value of 3 proposed by Wu through the theoretical analysis of the Korteweg-de Vries (KdV) equations. Finally, the phenomena of wave splitting and nonlinear focusing of a solitary wave propagating over a three-dimensional semicircular shoal are simulated. The results obtained from both the FNWD and WNWD models showing the fission process of separating a main solitary wave into multiple waves of decreasing amplitudes are presented, compared, and discussed.