The coupled motion-between multiple inviscid, incompressible, immiscible fluid layers in a rectangular vessel with a rigid lid and the vessel dynamics-is considered. The fluid layers are assumed to be thin and the shallowwater assumption is applied. The governing form of the Lagrangian functional in the Lagrangian particle path (LPP) framework is derived for an arbitrary number of layers, while the corresponding Hamiltonian is explicitly derived in the case of two-and three-layer fluids. The Hamiltonian formulation has nice properties for numerical simulations, and a fast, effective and symplectic numerical scheme is presented in the two-and three-layer cases, based upon the implicit-midpoint rule. Results of the simulations are compared with linear solutions and with the existing results of Alemi Ardakani et al. (J Fluid Struct 59:432-460, 2015) which were obtained using a finite volume approach in the Eulerian representation. The latter results are extended to non-Boussinesq regimes. The advantages and limitations of the LPP formulation and variational discretization are highlighted.