Abstract. In typical industrial practice based on IEC standards, wind turbine simulations are computed in the time domain for each mean wind speed bin using a few unsteady wind seeds. Software such as FAST, BLADED, or HAWC2 can be used to capture the unsteadiness and uncertainties of the wind in the simulations. The statistics of these aeroelastic simulation outputs are extracted and used to calculate fatigue and extreme loads on the wind turbine components. The minimum requirement of having six seeds does not guarantee an accurate estimation of the overall statistics. One solution might be running more seeds; however, this will increase the computation cost. Moreover, to move beyond blade element momentum (BEM)-based tools toward vortex/potential flow formulations, a reduction in the computational cost associated with the unsteady flow and uncertainty handling is required. This study illustrates the unsteady wind aerodynamic statistics' stationary character based on the standard turbulence models. This character is shown based on the output of National Renewable Energy Lab (NREL) 5MW reference machine BEM simulations. Afterwards, we propose a non-intrusive polynomial chaos expansion (PCE) to build a surrogate model of the loads' statistics, the rotor thrust, and torque, at each time step, to estimate the extreme statistics more accurately and efficiently.