Despite the increased availability of serological data, understanding serodynamics remains challenging. Serocatalytic models, which describe the rate of seroconversion (gain of antibodies) and seroreversion (loss of antibodies) within a population, have traditionally been fit to cross-sectional serological data to capture long-term transmission dynamics. However, a key limitation is their binary assumption on serological status, ignoring heterogeneity in optical density levels, antibody titers, and/or exposure history. Here, we implemented Gaussian mixture models - an established statistical tool - to cross-sectional data in order to characterize serological diversity of seasonal human coronaviruses (sHCoVs) across a wide range of age groups. These methods consistently identified multiple distinct seropositive levels, suggesting that among seropositive individuals, the number of prior exposures or response to infection may vary. We fit adapted, multi-compartment serocatalytic models with different assumptions on exposure history and waning of antibodies. The best fit model for each sHCoV was always one that accounted for host variation in the scale of serological response to infection. These models allowed us to estimate the strength and frequency of serological responses, finding that the time for a seronegative individual to become seropositive ranges from 2.40 to 7.03 years across sHCoVs, and most individuals mount a strong antibody response reflected in high optical density values, skipping lower levels of seropositivity. We find that despite frequent infection and strong serological responses, for all sHCoVs except 229E, individuals are likely to become seronegative again at some point after their first infection. Nonetheless, our results also indicate that by age 22, for each sHCoV the probability of having seroconverted at least once is over 95%. Crucially, our reimagined serocatalytic methods can be flexibly adapted across pathogens, having the potential to be broadly applied beyond this work.