Accurate estimation of extreme cold temperature return periods and their associated spatial uncertainties is critical for infrastructure design in cold climates. Traditional kriging interpolation provides optimal point predictions but fails to quantify spatial uncertainty, limiting risk-informed infrastructure design. We apply a five-step geostatistical workflow to quantify spatial uncertainty in extreme cold temperature ( \(T_{\min } \le -20, -30, -35, -40, -45\,^\circ \) C) return periods across Quebec’s 1.7 million \(\hbox {km}^2\) for the 1991–2020 climate normal period. The methodology integrates Gaussian transformations, linear regression revealing that latitude explains 28–46% of spatial variance, multivariate variogram modeling with elevation as a collocated covariate, multivariate Turning Bands Simulation generating 100 equiprobable realizations, and comprehensive uncertainty quantification through ensemble analysis. Results reveal pronounced geographical patterns in return period distributions. For moderate thresholds (− 30 °C), northern Quebec exhibits mean return periods of 2–3 years, while southern regions show return periods exceeding 20–30 years. For extreme thresholds (− 40 °C), events become rare even in the north (8–15 year return periods). Probability maps show a 70–90% likelihood of experiencing at least one − 40 °C event within any 10-year period in northern Quebec, compared to < 20% in southern areas. Uncertainty quantification reveals mean standard deviations ranging from 0.89 years (− 30 °C) to 10.04 years (− 40 °C), with maximum uncertainty in the central transition zone (48–52° N). These probabilistic maps enable risk-based infrastructure design. While this analysis characterizes the 1991–2020 baseline, observational and projected trends suggest a reduction in extreme cold frequency over coming decades, motivating extension of this framework to temporally non-stationary conditions.
