Cracks are common in nature and affect a number of soil properties. They increase the soil permeability, reduce the soil strength and stability, and thus affect the growth of plant roots. Characterization of crack networks in clay soils is needed to quantify the soil-plant relationships. This paper presents a numerical model of the desiccation cracking process, which is based on a lattice of Hookean springs with finite strength and critical strain. A new parameter, i.e., correlation length of critical strain, was used to reflect the spatially-correlated heterogeneity of the soil and to construct a critical strain random field. Three morphological parameters, which are indicators of quantitative crack analysis (i.e., area densities, length densities, and Euler number), were used to describe the spatial distribution and morphological connectivity of cracks. Field experimental results were used to evaluate the numerical model. The results showed that the model replicates the dynamic development of the crack network in nature and the general trends of shrinkage due to drying and crack development observed in the field. In addition, crack formation dynamics and final crack patterns were evaluated during the sensitivity analysis of correlation length. The model’s simulation of surface soil cracking was effective in predicting the formation of farmland soil cracks in clay soils.