Townsend’s model of attached eddies for boundary layers is revisited within a quasi-linear approximation. The velocity field is decomposed into a mean profile and fluctuations. While the mean is obtained from the nonlinear equations, the fluctuations are modelled by replacing the nonlinear self-interaction terms with an eddy-viscosity-based turbulent diffusion and stochastic forcing. Under this particular approximation, the resulting fluctuation equations remain linear, enabling solutions to be superposed, the same theoretical idea used in the original attached eddy model. By leveraging this feature, the stochastic forcing is determined self-consistently by solving an optimisation problem which minimises the difference between the Reynolds shear stresses from the mean and fluctuation equations, subject to a constraint that the averaged Reynolds shear-stress spectrum is sufficiently smooth in the spatial wavenumber space. The proposed quasi-linear approximation is subsequently applied to channel flow for Reynolds number Re𝜏 ranging from 500 to 20 000. The best result is obtained when the Reynolds stress is calculated by retaining only the two leading proper orthogonal decomposition modes, which further filters out the modelling artefact caused by the unphysical stochastic forcing. In this case, the resulting turbulence intensity profile and energy spectra exhibit the same qualitative behaviour as direct numerical simulation (DNS) data throughout the entire wall-normal domain, while reproducing the early theoretical predictions of the original attached eddy model within a controlled approximation to the Navier–Stokes equations. Finally, the proposed quasi-linear approximation reveals that the peak streamwise and spanwise turbulence intensities may deviate slightly from the logarithmic scaling with the Reynolds number for Re𝜏 ≳ 10 000, and the supporting evidence is presented using the existing DNS data.