Multiple scattering can severely affect the accuracy of optical instrumentation. Variance reduction methods have been implemented to improve a Monte Carlo model developed to simulate volume scattering functions measured by LISST-VSF instruments. The implemented methods can result in more than a tenfold increase in efficiency. The simulation is used to analyze multiple scattering errors for a range of Fournier-Forand (FF) phase functions. Our results demonstrate significant errors in the scattering coefficient, backscattering coefficient and phase function, where multiple scattering errors may only be considered negligible (<10%) for scattering coefficients <1 m-1. The errors depend strongly on the scattering coefficient but also increase when phase functions become more forward-peaked.