Catastrophic failures of partially or fully submerged structures, e.g., offshore platforms, hydrokinetic turbine blades, bridge decks, etc., due to the dynamic impact of free surface flows such as waves or floods have revealed the need to evaluate their reliability. In this respect, an accurate estimation of hydrodynamic forces and their relationship to instability in structures is required. The computational fluid dynamics (CFD) solver is known as a powerful tool to identify dynamic characteristics of flow; however, it commonly consumes a huge computational cost, especially in cases of re-simulations needed. In this paper, an efficient surrogate model based on the Gaussian process is developed to rapidly predict the nonlinear hydrodynamic pressure coefficients on submerged bodies near the water surface. For this purpose, a CFD model is first developed, which is based on a two-dimensional incompressible Navier–Stokes solver incorporating free surface treatment and turbulent flow models. Then, an experimental design is adopted to generate initial training samples considering the effect of the submerged body shape ratio and flow Re number. Surrogate models of hydrodynamic pressure coefficients and their instability based on Gaussian process modeling are established using the outcome from the CFD simulations, where optimal trend and correlation functions are also investigated. Once surrogate models are obtained, the mean and oscillation amplitudes of hydrodynamic pressure coefficients on a submerged rectangular body, which represents the shape of most civil structures, can be rapidly predicted without the attempt at re-simulation. The findings can be practically applied in rapidly assessing hydrodynamic forces and their instability of existing submerged civil structures or in designing new structures, where a suitable shape ratio should be adopted to avoid flow-induced instability of hydrodynamic forces.