Desiccation cracking, a common natural phenomenon, significantly affects the mechanical and hydraulic properties of clayey soils. This study employs 3D Discrete Element Method (DEM) to investigate the initiation and propagation mechanisms of soil desiccation cracking. By integrating water evaporation at the soil surface and water transportation between soil grains into a traditional DEM model, we unveil a sequential and hierarchical pattern in the development of desiccation cracks. Initially, micro cracks proliferate and coalesce in regions of soil defects characterized by the lowest coordination numbers, leading to the genesis of primary cracks. Subsequently, a complex network of cracks propagates, governed by the interplay between cracking and the redistribution of the tensile stress field. As desiccation advances, the detachment of soil from the bottom wall and an increase in soil strength contribute to the stabilization of the crack network. Notably, the hierarchical nature of desiccation cracking is predominantly shaped by initial soil defects and subsequently refined by stress concentration within the soil clods, correlating with changes in the length-to-height ratio of the soil clods.