This paper presents two methods for estimating the effect of numerical discretization error on statistical output accuracy in chaotic unsteady flow simulations: (1) an extension of recent advances in least-squares shadowing sensitivity calculations (LSS), and (2) an approximate time-windowing approach with individual adjoint solutions on each time window. Both methods rely on output adjoints, a direct application of which is not possible for chaotic systems, in which sensitivities to initial conditions and single-point discretization errors grow exponentially. This paper shows results for two prototypical chaotic systems: the Lorenz oscillator and the modified ergodic Kuramoto-Sivashinksy partial differential equation (MEKS). In addition it presents results for a low-Reynolds number Navier-Stokes flow to demonstrate the effectivity of the error estimates. Preliminary adaptive results are also included, in which spatially-localized forms of the error estimates drive static mesh adaptation to reduce errors in statistical outputs.