The eddy-covariance (EC) technique is considered the most direct and reliable method to calculate flux exchanges of the main greenhouse gases over natural ecosystems and agricultural fields. The resulting measurements are extremely important to characterize ecosystem exchanges of carbon, water, energy and other trace gases, and are widely used to validate or constrain parameter of land surface models via data assimilation techniques. For this purpose, the availability of both complete half-hourly flux time series and its associated uncertainty is mandatory. However, uncertainty estimation for EC data is challenging because the standard procedures based on repeated sampling are not suitable for this kind of measurements, and the presence of missing data makes it difficult to build any sensible time series model with time-varying second-order moments that can provide estimates of total random uncertainty. To overcome such limitations, this paper describes a new method in the context of the strategy based on the model residual approach proposed by Richardson et al. (Agric For Meteorol 148(1): 38–50, 2008). The proposed approach consists in (1) estimating the conditional mean process as representative of the true signal underlying observed data and (2) estimating the conditional variance process as representative of the total random uncertainty affecting EC data. The conditional mean process is estimated through the multiple imputation algorithm recently proposed by Vitale et al. (J Environ Inform 10.3808/jei.201800391, 2018). The conditional variance process is estimated through the stochastic volatility model introduced by Beltratti and Morana (Econ Notes 30(2): 205–234, 2001). This strategy is applied to ten sites that are part of FLUXNET2015 dataset, selected in such a way to cover various ecosystem types under different climatic regimes around the world. The estimated uncertainty is compared with estimates by other well-established methods, and it is demonstrated that the scaling relationship between uncertainty and flux magnitude is preserved. Additionally, the proposed strategy allows obtaining a complete half-hourly time series of uncertainty estimates, which are expected to be useful for many users of EC flux data.