Comparisons of genomes from recently diverged butterfly populations along a suture zone in central Texas have revealed high levels of divergence on the Z chromosome relative to autosomes, as measured by fixation index, $F_{st}$. The pattern of divergence appears to result from accumulation of incompatible alleles, obstructing introgression on the Z chromosome in hybrids. However, it is unknown whether this mechanism is sufficient to explain the data. Here, we simulate the effects of hybrid incompatibility on interbreeding butterfly populations using a model in which populations accumulate cross–incompatible alleles in allopatry prior to contact. We compute statistics for introgression and population divergence during contact between model butterfly populations and compare them to statistics obtained for 15 pairs of butterfly species interbreeding along the Texas suture zone. For populations that have evolved sufficiently in allopatry, the model exhibits high levels of divergence on the Z chromosome relative to autosomes in populations interbreeding on time scales comparable to periods of interglacial contact between butterfly populations in central Texas.Levels of divergence on the Z chromosome increase when interacting groups of genes are closely linked, consistent with interacting clusters of functionally related genes in butterfly genomes. Results for various periods in allopatry are in qualitative agreement with the pattern of data for butterflies, supporting a picture of speciation in which populations are subjected to cycles of divergence in glacial isolation, and partial fusion during interglacial contact.