This study proposes a new approach that can be used to generate the optimal surface state information and associated uncertainties from the estimates provided by the six land surface models used by the Global Land Data Assimilation System (GLDAS). The Förstner and best quadratic unbiased variance component estimators are used simultaneously with the least-squares method to calculate optimal values and the associated uncertainties. To demonstrate the concept, the research focused on three GLDAS hydrological products, namely soil moisture (SM), snow water equivalent (SWE), and canopy water (CAN) over the Canadian Prairies. When the Förstner estimator is applied, the estimated SM and SWE differ from their corresponding mean values by 26 mm and 9 mm respectively. Almost similar result was found with the best quadratic estimator. The estimated maximum uncertainties of each component including SM, SWE and CAN vary from year to year (e.g. 35 mm in 2006, 12 mm in 2007 and 2009 and 0.1 mm in 2001, respectively). The uncertainties of the total water storage (TWS) are almost similar to that of SM, which contributes more importantly to TWS in the area considered. The results obtained by the two proposed estimators were compared to the waterGAP hydrological models (WGHM), and to the Gravity Recovery and Climate Experiment (GRACE) terrestrial water storage anomalies. The optimal SWE anomalies generated from GLDAS using the proposed approach show a maximum correlation of r = 0.97 with the WGHM SWE anomalies. The optimal TWS anomalies have a correlation of r = 0.91 with WGHM, and r = 0.71 with GRACE. However, the correlation jumps to r = 0.81 when GRACE TWS is corrected for groundwater signals (with a mean RMSE of 8.5 mm). The RMSE and mean absolute error between our proposed methods and WGHM and GRACE are better than those obtained with each individual LSM or their average value. No significant mean bias error is observed in each case. Finally, the analysis of the time-lag characteristics of the resonance period between the results and their coherence was done by using a cross wavelet transform and a wavelet coherence analysis.