The multi-medium fluid flow around a supersonic sea-skimming flight is featured by the detached/attached shock wave, separation shock wave, and the reflected wave from the free surface. The complex wave structure and high nonlinearity pose a great challenge in accurate and stable numerical simulation. In this paper, a numerical model based on the high-order Runge–Kutta discontinuous Galerkin method is established to resolve the above problem. Based on the fact that the dimensionless vertical velocity is small, the air–water interface is linearized and a modified flux scheme is proposed to simplify the treatment of the multi-medium problem. A block-based adaptive mesh refinement scheme is adopted to capture the complex wave structure with the new nodes projected on the curved boundary. Finally, the numerical simulation of supersonic sea-skimming flight of the National Advisory Committee for Aeronautics 0012 airfoil is carried out by using the above-mentioned simplified numerical model based on the scheme of partition solution. The results show that the model can perform high-resolution simulations for the shock wave structure in various scenes. Meanwhile, the Mach number and distance between the airfoil and free surface are important factors affecting the structural characteristics of the shock wave systems and the airfoil loading characteristics. When the reflected shock wave acts on the airfoil's lower boundary, there will be a positive moment effect to make the airfoil dive, and the occurrence of this dangerous scene should be avoided. The relevant conclusions obtained can provide a reference for further research and engineering design.