High-frequency asymptotics is now a standard tool for analyzing localized hydrodynamic instabilities in many physical situations: flows with high Reynolds number, stratification, and rotation, to name a few. In general, the method gives an approximation which is first-order accurate in the asymptotic parameter. Here, a second-order accurate WKBJ (Wentzel–Kramers–Brillouin–Jeffreys) approximation is derived and its numerical properties are reported when applied to incompressible viscous flows. Numerical experiments compare the first- and second-order WKBJ approximations with direct numerical simulation (DNS) when computing the evolution of wave packet disturbances superimposed on a base flow. We analyze WKBJ’s convergence rate, time transient properties, and long-time behavior for eigenvalue calculations. Numerical experiments are performed with two base flows: Taylor–Green vortex, which is symmetric, and another vortex where the symmetry is broken. We find that the second-order WKBJ is excellent for approximating the transient features of wave packets in both base flows including the peaks of velocity related to the Orr mechanism of transient growth. Despite its success with transients, the second-order WKBJ remains accurate until some instant of time when, usually, it diverges from the DNS solution. Scale separation, essential for WKBJ, breaks down at this moment. We have identified a particular feature in the Lagrangian map of particles that explains, and predicts, where and when WKBJ and DNS diverge significantly. Because WKBJ may become imprecise when applied for long-time integration of viscous flows, it may yield incorrect approximation for eigenvalue calculations. Understanding its limitations is the key for successful application of WKBJ approximation.