A multi-layer hydrostatic shallow-water model was developed in the present study. The layer-integrated hydrostatic nonlinear shallow-water was solved with θ time integration and the least-squares finite element method. Since the least-squares formulation was employed, the resulting system of equations was symmetric and positive–definite; therefore, it could be solved efficiently by the preconditioned conjugate gradient method. The model was first applied to simulate the von Karman vortex shedding. A well-organized von Karman vortex street was reproduced. The model was then applied to simulate the Kuroshio current-induced Green Island vortex street. A swirling recirculation was formed and followed by several pairs of alternating counter-rotating vortices. The size of the recirculation, as well as the temporal and spatial scale of the vortex shedding, were found to be consistent with ADCP-CDT measurements, X-band radar measurements, and analysis of the satellite images. It was also revealed that Green Island vortices were affected by the upstream Orchid Island vortices.