A probabilistic model is presented to predict the formation and propagation of crack networks in thermal fatigue. It is based on a random distribution of sites where cracks initiate and on the shielding phenomenon corresponding to the relaxed stress field created around cracks. Currently, the model considers only heterogeneous uniaxial stress loading even if thermal fatigue is multiaxial. However, the first simulations on a uniaxial mechanical loading representative of the stress gradient that appears in thermal fatigue shocks are in qualitative agreement with experimental results. The larger the stress amplitude, the denser the crack network and the smaller the crack sizes.