Adaptive divergence and speciation may happen despite opposition by gene flow. Identifying the genomic basis underlying divergence with gene flow is a major task in evolutionary genomics. Most approaches (e.g. outlier scans) focus on genomic regions of high differentiation. However, not all genomic architectures potentially underlying divergence are expected to show extreme differentiation. Here, we develop an approach that combines hybrid zone analysis (i.e. focuses on spatial patterns of allele frequency change) with system-specific simulations to identify loci inconsistent with neutral evolution. We apply this to a genome-wide SNP set from an ideally-suited study organism, the intertidal snail Littorina saxatilis, which shows primary divergence between ecotypes associated with different shore habitats. We detect many SNPs with clinal patterns, most of which are consistent with neutrality. Among non-neutral SNPs, most are located within three large putative inversions differentiating ecotypes. Many non-neutral SNPs show relatively low levels of differentiation. We discuss potential reasons for this pattern, including loose linkage to selected variants, polygenic adaptation and a component of balancing selection within populations (which may be expected for inversions). Our work is in line with theory predicting a role for inversions in divergence, and emphasises that genomic regions contributing to divergence may not always be accessible with methods purely based on allele frequency differences. These conclusions call for approaches that take spatial patterns of allele frequency change into account in other systems.