The main problem about a rigorous solution to such a scattering proplem is that computations are extremely demanding. Just imagine you have a wavelength $\lambda$ of some $400$nm to $700$nm for visible light (from here):
Now, to do physically meaningful simulations, you will need a sub-wavelength lattice which makes any computational cell above, say $10\,\mu m^3$ not accessible since you have in the order of one million grid points.
Approximative Approaches
But of course there can be ways out of it if you are willing to make some approximations which will largely depend on the characteristics of the particles you are looking at. It is best to assume that we only have spherical particles since we can apply Mie theory in this case.
Large Particles
First of all, let us consider particles which are much larger than the wavelength. Then, the radius $R$ times the wave vector $k=2\pi/\lambda$ is much bigger than one, $$kR\gg1$$
which basically means that one observes reflection at a plane interface. You can implement these particles using geometrical optics (mixed with Fresnel reflection if you like) since nothing really wave-like will happen as in this image (taken from here):
Small Particles
Second, the particles should be much smaller than the wavelength,
$$kR\ll1\,.$$
Then, everything what is observed is a sum of dipolar responses of the particles in the so-called Rayleigh-scattering. Then,
the intensity of light scattered by a single small particle from a beam of unpolarized light of wavelength $\lambda$ and intensity $I_0$ is given by:
$$I=I_0(1+\cos^2\theta)\frac{(kR)^6}{2(kr)^2}\left(\frac{n_p^2-1}{n_p^2+2}\right)$$
where I have chosen the variables to be consistent with the used terminology and $r$ is the distance to the object, $\theta$ is the scattering angle and $n_p$ is the sphere's refractive index. Here is an image of such a situation with some metal particles also having quadrupolar excitation (from here):
A Mean Field Approach - Effective Permittivity
If you have a lot of these small objects, you may use the Clausius-Mossotti relation which gives you an effective permittivity $\epsilon_p=n_p^2$ depending on the concentration of the particle in some volume:
$$\epsilon_{eff} = \epsilon_p + \frac{n\alpha}{1-\frac{n\alpha}{3\epsilon_p}}$$
where $\alpha$ is the polarizability of the sphere, for details see e.g. Electromagnetic mixing formulas and applications by Sihvola. This would be something like a mean-field approach. You can make some very neat effects using this effective approach since it allows you to calculate a continuous refraction around some particle streams under water.
However, if the particles size is in the order of the wavelength, $$kR\approx 1$$ then you may have to take higher multipole moments into account which may be a very demanding task.
For much more on the subject I would recommend Bohren & Huffmanns classic Absorption and Scattering of Light by Small Particles.
Sincerely
