Tsunamis entering into shallow water regions may become highly nonlinear and this may be due to the irregularity of sea bottom roughness relative to the water depth and the complex coastline geometry. The elliptic mild-slope equation is commonly used to predict the nonlinear wave propagation in shallow water regions but it requires huge amount of computer resources which may not be practical for tsunami propagation predictions. An efficient finite element approach has been adopted in the present project to resolve the nonlinear problem of wave transformation in near shore zones as well as to better conform the model grids to any complex coastline configurations. The efficient approach is based on the wave action conservation equation that takes into account of wave refraction-diffraction and energy dissipation due to bottom roughness. An operator splitting scheme is employed to solve the wave action conservation equation. Firstly, to increase numerical stability, the Eulerian-Lagrangian method is applied to solve the advection terms in the equation. The horizontal terms are then discretized by an implicit finite element method and, finally, the vertical terms are approximated by an implicit finite difference method. A nominal-time iteration method is used to efficiently solve the non-linear irrotational wave number equation for the wave direction. Over 6000 nine-node elements have been used to mesh the Pearl River estuary region. The boundary conditions are based on the results obtained from a simulation applied for a larger computation domain encompassing the entire South China Sea. The computed result provides a general picture of tsunami propagation in the desired region. Model validation and result verification, however, are necessary for any future prediction exercises.