This article describes development of a model to predict microcrack density evolution as a function of thermal history in porous ceramics with anisotropic thermal expansion coefficients. The model is based on the evolution of statistical distributions of grain-boundary stresses and crack opening displacements (CODs), and represents universal cracking and crack closure behavior scaled by coefficients that depend on the microstructure and material properties. The functional forms for this model were developed from a 2D finite element simulation of thermally induced microcracking. Specifically, the grain-boundary stresses and CODs were studied, and a parameter describing each quantity was calibrated from modulus hysteresis data. An additional calibration parameter was also introduced to link crack density to macroscopic behavior. The model was applied to experimental results from a porous cordierite ceramic and shown to accurately describe its behavior.