In this paper, we will combine an upwind radial basis function-finite element with direct velocity–pressure formulation to study the two-dimensional Navier-Stokes equations with free surface flows. We will examine this formulation in an improved mixed-order finite element and localized radial basis function method. A particle tracking method and the arbitrary Lagrangian-Eulerian scheme will then be applied to simulate the two-dimensional high Reynolds free surface flows. An upwind improved finite element formulation based on a localized radial basis function differential quadrature (LRBFDQ) method is used to deal with high Reynolds number convection dominated flows. This study successfully obtained very high Reynolds number free surface flows, up to Re = 500 000. Finally, we will demonstrate and discuss the capability and feasibility of the proposed model by simulating two complex free surface flow problems: (1) a highly nonlinear free oscillation flow and (2) a large amplitude sloshing problem. Using even very coarse grids in all computing scenarios, we have achieved good results in accuracy and efficiency.