The boundary element method (BEM) has been utilized to formulate a three-dimensional math model for the foam earplug-earcanal and earmuff-earcanal systems. The BEM technique consists of transforming the partial differential equation describing the behavior of the unknowns, inside and on the boundary of the domain, into an integral equation relating only boundary values. In representing the viscoelastic properties of the foam earplug, both integer differential operator and fractional operator constitutive equations were utilized. The BEM model was then utilized to study the steady-state and transient responses (including the prediction of each system's insertion loss) for different earplug-earcanal and earmuff-earcanal configurations. The basic BEM model has demonstrated the ability to predict the internal resonances for a production earmuff shell-earcanal system and the measured insertion loss of a foam earplug-earcanal system subjected to an impulse. Together with these findings a general discussion of the need for improvements in the ability to model a variety of earplug- and earmuff-earcanal configurations will be presented.