We propose an efficient multi-level Monte Carlo ensemble domain decomposition method for solving random Stokes-Darcy coupled models with uncertain hydraulic conductivity. The multi-level Monte Carlo (MLMC) method is employed to reduce the computational cost in the probability space by decreasing the number of simulations with spatial mesh refinement. Meanwhile, an ensemble domain decomposition approach is introduced to improve efficiency in the physical space, allowing multiple samples to share a common coefficient matrix.

We develop a finite element framework for the proposed method and analyze its convergence properties under different choices of Robin interface parameters. In particular, we establish both mesh-dependent and mesh-independent convergence results and provide optimized parameter selections to accelerate the iterative procedure. Numerical experiments are presented to validate the theoretical analysis and demonstrate the efficiency of the proposed method.