A probabilistic model is proposed to predict the formation and propagation of crack networks in thermal fatigue. It is based on a random distribution of sites where cracks form and on the shielding phenomenon corresponding to the relaxed stress field created around propagating cracks. The stress gradient that arises under a surface submitted to thermal shocks is accounted for as well as multiaxial stress states. Experiments using digital image correlation have been performed to introduce the hypotheses made herein, and others to identify the parameters of the crack density as a function of the stress range and the number of cycles. Simulations of the formation of crack networks in a heterogeneous biaxial state of stress are carried out. A good qualitative agreement is obtained with existing thermal fatigue experiments especially considering the effect of the stress range on the crack densities as well as on the distribution of crack sizes.