© 2015 Society for Industrial and Applied Mathematics. We present a mathematical model to describe the distribution of surfactant pairs in a multilayer structure beneath an adsorbed monolayer. A mesoscopic model comprising a set of ordinary differential equations that couple the rearrangement of surfactant within the multilayer to the surface adsorption kinetics is first derived. This model is then extended to the macroscopic scale by taking the continuum limit that exploits the typically large number of surfactant layers, which results in a novel third-order partial differential equation. The model is generalized to allow for the presence of two adsorbing boundaries, which results in an implicit free-boundary problem. The system predicts physically observed features in multilayer systems such as the initial formation of smaller lamellar structures and the typical number of layers that form in equilibrium.