Missing data is a ubiquitous problem in evaluating long-term experimental measurements, such as those associated with the FluxNet project, due to the equipment failures, system maintenance, power-failure, and lightning strikes among other things. To estimate annual values of net ecosystem carbon exchange (NEE), latent heat flux (LE) and sensible heat flux (H), such gaps in the measured data must be filled or imputed. So far, no standardized method has been accepted and the imputation methods used are largely dependent on the researchers’ choice. Here, we used multiple imputation (MI) to gap-fill the missing data for annual estimations of NEE, LE and H at three flux sites associated with the FluxNet effort. MI is a Monte Carlo technique in which the missing values are replaced by several simulated values. Each data set imputed is a complete one where the observed values are the same as those in the original data set; only the missing values are different. Thus, the normal statistical analysis (e.g. annual total calculation) can be applied to each data set separately. The results of each analysis can be recombined into one summary. We applied the MI method to eddy covariance measurements collected from Walker Branch Watershed (WBW) site (a deciduous forest), Duke site (a coniferous forest) and Niwot site (a subalpine forest). Results showed that annual estimations of NEE, LE and H by MI were comparable to other imputation methods but MI was much easier to apply because of readily available software and standard algorithms. Besides the normal statistical analyses, MI also provided confidence intervals for each estimated parameter. This confidence interval is most useful when assessing energy, water, and carbon balance closures at a given tower site. Significant differences in annual NEE, LE and H were found among years at the three AmeriFlux sites. NEE at the Niwot Ridge site was lower and LE and H were higher than at the other two sites. With the available software and realistic gap-filling capability, MI has the potential to become a standardized method to gap-fill eddy covariance flux data for annual estimations and to improve the analysis of uncertainties associated with annual estimations of NEE, LE and H from regional and global flux networks.