A three-dimensional numerical model for solving JONSWAP waves and jets in σ-coordinate was developed using large eddy simulation technology, in which the dynamic coherent eddy model was used. The operator splitting method was used to solve the Navier-Stokes equations. After the model was verified with the experimental results of the jet in regular waves, the model was used to simulate the vertical turbulent jet in JONSWAP wave environments. Investigation was done for the distribution of the enstophy of the core of vorticity in the backzone. The mixing process was found to be governed by entrainment and jet-wave interaction. The mechanism of jet instantaneous characteristic in JONSWAP wave was also investigated. The numerical results are of great theoretical importance for the understanding of jet turbulent behaviors in JONSWAP wave.