A numerical technique, based on the boundary integral method, is developed to allow the modelling of unsteady free-surface flows at large Reynolds numbers in cases where the surface is contaminated by some surface-active compound. This requires the method to take account of the tangential stress condition at the interface and is achieved through a boundary-layer analysis. The constitutive relation that forms the surface stress condition is assumed to be of the Boussinesq type and allows the incorporation of surface shear and dilatational viscous forces as well as Marangoni effects due to gradients in surface tension. Sorption kinetics can be included in the model, allowing calculations for both soluble and insolube surfactants. Application of the numerical model to the problem of bursting gas bubbles at a free surface shows the greatest effect to be due to surface dilatational viscosity which drastically reduces the amount of surface compression and can slow and even prevent the information of a liquid jet. Surface tension gradients give dilatational elasticity to the surface and thus also significantly prevent surface compression. Surface shear viscosity has a smaller effect on the interface motion but results in initially increased surface concentrations due to the sweeping up of surface particles ahead of the inward-moving surface wave.