The objective of this paper is to present an efficient beam formulation based on the boundary element method (BEM), for the nonlinear dynamic analysis of wind turbine towers of variable cross section founded either on surface or on monopile foundation system. The whole structure may undergo moderately large displacements, taking into account the effect of soil–structure kinematic and inertia interaction. The tower is subjected to the combined action of arbitrarily distributed or concentrated transverse wind loading as well as to seismic excitation together with axial loading arising from the self-weight of the tower and the mechanical parts. The Blade element momentum theory is taken into consideration in order to produce the wind load time histories, while the site seismic response is obtained through one dimensional shear wave propagation analysis. The soil–surface foundation system is formulated as equivalent lateral and rotational springs, while the case of monopile system is treated as a prismatic beam on elastic foundation assigning the corresponding springs and dashpots along its length. An extensive case study is carried out on a wind turbine tower–foundation system employing the generated wind velocity time histories and recorded earthquake accelerograms, providing insight to several structural phenomena. The results of the proposed model are compared wherever possible with those obtained from a commercial finite elements software package, illustrating the validity and the efficiency of the developed method. From the obtained results, the strong influence of the nonlinear effects on the dynamic response of the wind turbine tower is verified.