The joint likelihood of observable cluster signals reflects the astrophysical evolution of the coupled baryonic and dark matter components in massive halos, and its knowledge will enhance cosmological parameter constraints in the coming era of large, multi-wavelength cluster surveys. We present a computational study of intrinsic covariance in cluster properties using halo populations derived from Millennium Gas Simulations (MGS). The MGS are re-simulations of the original 500 Mpc/h Millennium Simulation performed with gas dynamics under two different physical treatments: shock heating driven by gravity only (GO) and a second treatment with cooling and preheating (PH). We examine relationships among structural properties and observable X-ray and Sunyaev-Zel'dovich (SZ) signals for samples of thousands of halos with M_200 > 5 \times 10^{13} Msun/h and z < 2. While the X-ray scaling behavior of PH model halos at low-redshift offers a good match to local clusters, the model exhibits non-standard features testable with larger surveys, including weakly running slopes in hot gas observable--mass relations and ~10% departures from self-similar redshift evolution for 10^14 Msun/h halos at redshift z ~ 1. We find that the form of the joint likelihood of signal pairs is generally well-described by a multivariate, log-normal distribution, especially in the PH case which exhibits less halo substructure than the GO model. At fixed mass and epoch, joint deviations of signal pairs display mainly positive correlations, especially the thermal SZ effect paired with either hot gas fraction (r=0.88/0.69 for PH/GO at z=0) or X-ray temperature (r=0.62/0.83). We discuss halo mass selection by signal pairs, and find a minimum mass scatter of 4% in the \PH model by combining thermal SZ and gas fraction measurements.