  We analyze a fully discrete finite element numerical scheme for the
Cahn-Hilliard-Stokes-Darcy system that models two-phase flows in coupled free
flow and porous media. To avoid a well-known difficulty associated with the
coupling between the Cahn-Hilliard equation and the fluid motion, we make use
of the operator-splitting in the numerical scheme, so that these two solvers
are decoupled, which in turn would greatly improve the computational
efficiency. The unique solvability and the energy stability have been proved
in~\cite{CHW2017}. In this work, we carry out a detailed convergence analysis
and error estimate for the fully discrete finite element scheme, so that the
optimal rate convergence order is established in the energy norm, i.e.,, in the
$\ell^\infty (0, T; H^1) \cap \ell^2 (0, T; H^2)$ norm for the phase variables,
as well as in the $\ell^\infty (0, T; H^1) \cap \ell^2 (0, T; H^2)$ norm for
the velocity variable. Such an energy norm error estimate leads to a
cancellation of a nonlinear error term associated with the convection part,
which turns out to be a key step to pass through the analysis. In addition, a
discrete $\ell^2 (0;T; H^3)$ bound of the numerical solution for the phase
variables plays an important role in the error estimate, which is accomplished
via a discrete version of Gagliardo-Nirenberg inequality in the finite element
setting.
