Land surfaces dissipate energy through latent (LE) and sensible (H) heat fluxes that modulate atmospheric temperature and humidity, which in return affect land surface vegetation and soil processes. Within this two-way land-atmosphere coupling, surface energy partitioning (LE versus H) plays a central role in connecting the land and atmosphere states and fluxes. However, considerably large uncertainties still exist in earth system land models, i.e. the phase 6 of the Coupled Model Intercomparison Project (CMIP6). Further, the underlying controls from climate and biological factors on surface energy partitioning over different biome types are not well understood. In this study, we combined machine learning (ML) and causal inference models to investigate and reduce the uncertainties (i.e., parametric, structural, and forcing uncertainties) of CMIP6 simulated evaporative fraction (defined as LE / (LE + H)) across 64 FLUXNET sites that cover five major biomes. We found that CMIP6 model ensembles overestimated evaporative fraction with considerable spread at deciduous broadleaf forest, evergreen needleleaf forest, and savanna sites. Accounting for the biases from all related surface climate and biological driving variables, the CMIP6 model simulated EF could be largely improved (e.g., R between model and observed EF improved from 0.47 to 0.66), with leaf area index, vapor pressure deficit, and precipitation dominated the model improvement. Furthermore, ML-based parameterization generally showed a promising opportunity to further reduce model biases (e.g., R improved from 0.66 to 0.80) in spite of the limited improvement at evergreen broadleaf forest sites where model bias may be dominated by structural imperfection. This study provided an effective framework to understand and reduce model uncertainties in simulating land surface energy flux partitioning and, more importantly, highlighted the need of effective model structure improvement for the next generation earth system land model development.