Flow systems with highly nonlinear free/moving surface motion are common in engineering applications, such as wave impact and fluid-structure interaction (FSI) problems. In order to reveal the dynamics of such flows, as well as provide a reduced-order modeling (ROM) for large-scale applications, we propose a proper orthogonal decomposition (POD) technique that couples the velocity flow field and the level-set function field, as well as a proper normalization for the snapshots data so that the low-dimensional components of the flow can be retrieved with a priori knowledge of equal distribution of the total variance between velocity and level-set function data. Through numerical examples of a sloshing problem and a water entry problem, we show that the low-dimensional components obtained provide an efficient and accurate approximation of the flow field. Moreover, we show that the velocity contour and orbits projected on the space of the reduced basis greatly facilitate understanding of the intrinsic dynamics of the flow systems.