The problem of crack pattern formation due to thermal shock loading at the surface of half space is solved numerically using the two-dimensional boundary element method. The results of numerical simulations with 100–200 random simultaneously growing and interacting cracks are used to obtain scaling relations for crack length and spacing. The numerical results predict that such a process of pattern formation with quasistatic crack growth is not stable and at some point the excess energy leads to unstable propagation of one of the longest cracks. This single-crack scenario should be understood in a local sense. There could be other unstable cracks far away that together can form a new pattern. The onset of instability has also been determined from numerical results.