Numerical approximation of stochastic Stokes–Darcy equations usually requires repeated sampling of the random hydraulic conductivity tensor and then simulating flow ensembles. In this setting, we propose an efficient, second order, ensemble algorithm for fast computation of the whole set of realizations of the stochastic Stokes–Darcy model corresponding to different random hydraulic conductivity tensor samples. The ensemble algorithm only requires the solution of two linear systems that have the same constant coefficient matrices for all realizations. We give a complete long time stability and convergence analysis for the method. Numerical experiments are presented to support theoretical results and demonstrate the application of the method.