Damage to components made of brittle material due to thermal shock represents a high safety risk. Predicting the degree of damage is therefore very important to avoid catastrophic failure. An energy-based linear elastic fracture mechanics bifurcation analysis using a three-dimensional finite element model is presented here, which allows the determination of crack length and crack spacing for a defined thermal load in a free plate. It is assumed that a hierarchical crack pattern is formed due to cooling penetration. The constant growth of the ideal regular pattern of hexagons can change into a pattern with a different symmetry by slightly changing the cooling conditions. This bifurcation point is determined by the second derivative of the mechanical potential with respect to the geometry of the crack front. The very high computational effort for the second derivative is reduced by describing the three-dimensional crack front with a limited number of Fourier coefficients. A one-dimensional transient temperature field at a sufficient distance from the plate edge is assumed. For alumina, the crack length and crack spacing curves are computed for different quenching temperatures and heat transfer coefficients. The corresponding final crack lengths are also calculated as a measure of damage. Comparison with a two-dimensional model confirms the expected 1/2 difference in crack spacing. Data from thermal shock experiments are also presented. However, due to the cracks caused by the strong cooling at the edge, these correspond to the results of the two-dimensional model.