Enhancement of the FLARE method frequently employed in a parabolic-equation model analysis for buoyancy-induced flow reversal within vertical channels is investigated. The present study shows several modifications to the parabolic model and the FLARE method. The relative performance of all these modified parabolic models is evaluated for two-dimensional buoyancy-assisted reversed flows. Comparison of the solutions between the parabolic and the full elliptic Navier-Stokes models is also made. Results tend to confirm previously reported findings that analysis based on the parabolic model and the FLARE approximation ( B.L. + FLARE) reduces the implementation and computation efforts remarkably yet still retains the accuracy of the numerical solutions. However, when the Richardson number is higher than a threshold number, the traditional BX. + FLARE procedure leads to an oscillating or even a divergent solution downstream. The numerical instability downstream can be eliminated by introducing a modified numerical algorithm ( BX. + FLAREC method) proposed in this study so as to advance the analysis into the higher Richardson number regime. For example, considering the channel flow at Re = 500 and Ri = 0.9, the BX. + FLARE method leads to an oscillating solution in flow and temperature fields downstream; however, BX. + FLAREC can still provide a stable solution that closely agrees with the elliptic-model results.