Direct numerical simulations are conducted for temporally evolving stratified wake flows at Reynolds numbers from
$10\,000$
to
$50\,000$
and Froude numbers from
$2$
to 50. Unlike previous studies that obtained statistics from a single realization, we take ensemble averages among 80–100 realizations. Our analysis shows that data from one realization incur large convergence errors. These errors reduce quickly as the number of statistical samples increases, with the benefit of ensemble average diminishing beyond 40–60 realizations. The data with ensemble average allow us to test the previously established scalings and arrive at new scaling estimates. Specifically, the data do not support power-law scaling in the centreline velocity deficit
$U_0$
beyond the near wake. Its decay rate increases continuously from 0.1 at the onset of the non-equilibrium regime until the end of our calculations without reaching any asymptote. Additionally, while no power-law scalings could be found in the wake width (
$L_H$
) and wake height (
$L_V$
) in the late wake,
$L_H\sim (Nt)^{1/3}$
is a good working approximation of the wake's horizontal size, where
$N$
is buoyancy frequency and
$t$
is time. Besides the low-order statistics, we also report the transverse integrated terms and the vertically integrated terms in the turbulent kinetic energy budget equation as a function of the vertical and transverse coordinates. The data indicate that there are two peaks in the vertically integrated production and transport terms, and one peak when the two terms are integrated horizontally.