  We propose, analyze and test a fully discrete, efficient second-order
algorithm for computing flow ensembles average of viscous, incompressible, and
time-dependent magnetohydrodynamic (MHD) flows under uncertainties in initial
conditions. The scheme is decoupled and based on Els\"asser variable
formulation. The algorithm uses the breakthrough idea of Jiang and Layton, 2014
to approximate the ensemble average of $J$ realizations. That is, at each time
step, each of the $J$ realization shares the same coefficient matrix for
different right-hand side matrices. Thus, storage requirements and
computational time are reduced by building preconditioners once per time step
and reuse them. We prove stability and optimal convergence with respect to the
time step restriction. On some manufactured solutions, numerical experiments
are given to verify the predicted convergence rates of our analysis. Finally,
we test the scheme on a benchmark channel flow over a step and it performs
well.
