The paper proposed a two-phase flow Large Eddy Simulation (LES) algorithm for atomization in supersonic gas flow. A coupled Level Set and Volume of Fluid (VOF) interface tracking method is used. The supersonic gas flow is solved using a compressible flow solver, and the liquid phase is solved using an incompressible flow solver. Appropriate boundary conditions are specified on the interface for both solvers to correctly capture the interaction between the gas and liquid. The water jet primary breakup in a supersonic air crossflow is well predicted, with the shock waves ahead of the liquid column and drops/ligaments properly captured. Two kinds of vortices respectively relating to surface breakup and column breakup are observed. Surface waves on the windward side of the liquid jet are numerically reproduced, and the wavelength is well correlated to an effective Weber number which is determined by the gas flow behind the shock. As gas flow goes across the strong shock wave ahead of liquid column, the disintegrating aerodynamic force is weakened. Rayleigh-Taylor instability is found to be responsible for the development of surface waves which causes column breakup. It is demonstrated that LES can produce almost the same breakup morphology as an under-resolved DNS. (C) 2016 Elsevier Ltd. All rights reserved.