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üller 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. Only after a detailed numerical treatment of both diagonal and close-to-diagonal matrix elements is the unitarity condition found to be well satisfied for the nonabsorbing case (U>0.995), a result that testifies to the accuracy of our approach.