Nowadays, various fields in environmental sciences require the availability of appropriate techniques to exploit the information given by multivariate spatial or spatio-temporal observations. In particular, radon flux data which are of high interest to monitor greenhouse gas emissions and to assess human exposure to indoor radon are determined by the deposit of uranium and radio (precursor elements). Furthermore, they are also affected by various atmospheric variables, such as humidity, temperature, precipitation and evapotranspiration. To this aim, a significant role can be recognized to the tools of multivariate geostatistics which supports the modeling and prediction of variables under study. In this paper, the spatio-temporal distribution of radon flux densities over the Veneto Region (Italy) and its estimation at unsampled points in space and time are discussed. In particular, the spatio-temporal linear coregionalization model is identified on the basis of the joint diagonalization of the empirical covariance matrices evaluated at different spatio-temporal lags and is used to produce predicted radon flux maps for different months. Probability maps, that the radon flux density in the upcoming months is greater than three historical statistics, are then built. This might be of interest especially in summer months when the risk of radon exhalation is higher. Moreover, a comparison with respect to alternative models in the univariate and multivariate context is provided.