  An accurate and efficient numerical simulation approach to electromagnetic
wave scattering from two-dimensional, randomly rough, penetrable surfaces is
presented. The use of the M\"uller equations and an impedance boundary
condition for a two-dimensional rough surface yields a pair of coupled
two-dimensional integral equations for the sources on the surface in terms of
which the scattered field is expressed through the Franz formulas. By this
approach, we calculate the full angular intensity distribution of the scattered
field that is due to a finite incident beam of $p$-polarized light. We
specifically check the energy conservation (unitarity) of our simulations (for
the non-absorbing case). Only after a detailed numerical treatment of {\em
both} diagonal and close-to-diagonal matrix elements is the unitarity condition
found to be well-satisfied for the non-absorbing case (${\mathcal U}>0.995$), a
result that testifies to the accuracy of our approach.
