A probabilistic closure model to include the energy lost in surrogate models based on proper-orthogonal-decomposition (POD) modes is introduced. A transformer-encoder block is utilized to learn the spatial probability density function of fluctuations that were present in the truncated POD modes. A proof of concept is done on the wake of the Windsor body at yaw angles of δ = [2.5◦,5◦,7.5◦,10◦,12.5◦], with δ = 7.5◦ as a test case. The key coherent modes are identified by clustering them based on dominant frequency dynamics. To do so, Hotelling’s T2 is applied on the power spectral density of the POD temporal coefficients. These coherent modes account for nearly 60% of the total energy while comprising less than 10% of all modes and they would be the core of a surrogate model. A new POD basis is created by concatenating and orthonormalizing the coherent modes from training angles. This reduces the basis vectors from 142 to 90 without losing information. Transformers with different size on the attention layer, (64, 128 and 256), are trained to model the missing fluctuations. Larger attention sizes always improve predictions for the training set, but the transformer with an attention layer of size 256 overshoots the fluctuations predictions. Such overshoot is explained by the lower intensity of the fluctuations in the test set compared to the ones in the training cases. Adding the predicted fluctuations closes the energy gap between the reconstruction and the original flow field, improving predictions for energy, root-mean-square velocity fluctuations, and instantaneous flow fields. The deepest architecture reduces mean energy error from 37% to 12% and decreases the Kullback–Leibler divergence of velocity distributions from DK L = 0.2 to below DK L = 0.026.