Quantifying the structural uncertainty inherent in General Circulation Models (GCMs) is a fundamental prerequisite for robust climate change assessments. To mitigate these intrinsic uncertainties, which propagate errors throughout the simulation framework, employing multi-model ensembles (MMEs) is substantially more effective than relying on individual models. However, conventional arithmetic averaging within MMEs exhibits an inherent methodological limitation: by indiscriminately assigning equal statistical weight to every member, this unweighted approach obscures vital differences in model architecture and projected skill. Consequently, the composite output becomes highly susceptible to skewness driven by a minority of models exhibiting substantial systematic bias. The present study investigates 32 downscaled models from the NEX-GDDP-CMIP6 dataset over the historical period (1980–2014), benchmarking them against CRU TS observational records. Concurrently, we evaluate the capacity of Reliability Ensemble Averaging (REA) to constrain the uncertainty bounds of these modeled outputs. Evaluating the individual models, the unweighted ensemble, and the REA-derived construct against historical observations necessitates a multidimensional assessment framework. Isolated statistical metrics, such as basic correlation or mean bias, frequently fail to capture the complex error architecture inherent in GCM simulations. To bridge this gap, integrating Taylor diagrams with the non-parametric Kling-Gupta Efficiency (KGEnp) and the coefficient of determination (R2) provides the analytical depth needed to decouple systematic biases from purely random errors. Our analysis reveals that certain models severely overestimated surface air temperatures across both the high-altitude cold zones of the Zagros and Alborz mountain ranges and the hyper-arid central and southern plains, most notably over the Lut Desert and the southern coastal margins. Specific ensemble members, including CMCC-CM2-SR5, TaiESM1, and NorESM2-MM, generated considerable warm biases. Conversely, members of the GFDL, ACCESS, and MPI model families captured sharp spatial and topographic gradients better. While simple arithmetic averaging partially dampened these model-specific deviations, the unweighted MME still exhibited a persistent warm bias, overestimating the regional mean (17.6°C) relative to observations (17.4°C). The REA algorithm addresses this by assigning statistical weights at individual grid cells based on dual criteria: historical model performance and cross-model convergence. This localized weighting mechanism proved highly effective at suppressing systematic deviations. Specifically, the REA framework elevated R2 to 0.71 and KGEnp to 0.85, distinctly outperforming the simple unweighted ensemble (0.66 and 0.82, respectively). Beyond enhancing spatial correlation, this strategy appreciably reduced the uncertainty envelope and delineated high-altitude thermal boundaries with much sharper precision than standard averaging could achieve. These findings indicate that the REA algorithm provides a highly robust framework for accurately representing historical temperature profiles across Iran’s topographically complex and arid landscapes. Ultimately, filtering climate signals via performance-based weighting establishes a rigorous foundation for minimizing error variance. This underscores that the most effective way to reduce regional-scale uncertainty—and extract dependable climate signals for future projections—is through reliability-driven, weighted multi-model ensembles.