In two earlier papers, we studied the statistical and mechanistic structure of the turbulent boundary layer under a stress-free (clean) free surface. Findings there, such as the presence of inner and outer surface layers, are very much the direct result of the absence of shear stresses at the surface. The latter condition is easily lost when the surface is contaminated and surface elasticity varies with space and time. In this paper we consider the effect of surfactant on features of the free-surface turbulent flow. We perform direct numerical simulations of the Navier–Stokes equations subject to surfactant-laden free-surface boundary conditions for varying Reynolds and Marangoni numbers and low Froude numbers. As expected, the Marangoni effect decreases the horizontal turbulence intensity and normal vorticity at the surface. The direct effect on the turbulent kinetic energy is an increase in the dissipation and viscous diffusion and a decrease in the production near the surface relative to the clean case. The most prominent effect of the presence (of even a small amount) of surfactant is the drastic reduction in the surface divergence and the associated sharp decrease of up- and downwelling at the surface which has direct implications to near-surface turbulent transport. The observed surfactant effects on turbulent kinetic energy budget can be attributed to the generation of Marangoni vorticity at the free surface by approaching hairpin vortices. The Marangoni effect has also a direct effect on the boundary-layer structure, causing an increase of the thickness of the boundary layer and in the maxima of the mean shear near the surface. For moderate values of the Marangoni number, up-/downwelling effectively vanishes and the flow approaches a state independent of the Marangoni number. Guided by these results and to obtain theoretical insight, we develop a similarity solution for the mean flow. The analytic solution agrees well with the numerical data and provides precise measures for the multi-layer structure of the boundary layer. Based on the theoretical model, we derive scaling laws for the thickness of the inner and the outer boundary layers, which are also confirmed by numerical simulations.