Rayleigh scattering light intensity Assume a single ray of unpolarized light propagating through an opaque medium (water). Due to absorption and Rayleigh scattering the intensity decreases and light gets scattered/diffused before hiting the ground. Attenuation by absorption is given by
$$I = I_0 \exp(-\mu_{att}\cdot l_{depth}).$$
But, how does the light intensity profile on the ground, including scattering, look like? Formula?
Edit: Think the solution might be derived from this.
 A: Unlike attenuation, scattering is difficult because energy which is scattered out of the beam can be scattered back in.

If scattering coefficient is very weak over the path length of interest, or the scattering strongly favours small angular deviations, (green path in diagram) you can add the scattering and attenuation coefficients together to estimate the intensity at the end location, and use the Beer-Lambert law,
$$I_\nu(x)=I_{\nu,0} \exp\left(-\left(k_{\nu, s} + k_{\nu, a}\right)x\right)$$
Obviously, this only holds provided that light scattered out of the beam does not return. You will need to determine if this holds for your purpose, or not. 
In reality, there will always be in-scattering, and these effects can be strong (see red path in the diagram).
So the full answer is that you need to solve the radiative transfer equation. Assuming you can ignore any heating effects, and you just want to know about attention of the beam, you can set the emission coefficient to zero and find the solution of,
$$
\hat{\Omega} \cdot \nabla I_\nu + (k_{\nu, s}+k_{\nu, a}) I_\nu = \frac{1}{4\pi}k_{\nu, s} \int_\Omega I_\nu d\Omega
$$
However, there are whole books written on this topic. You will probably want to make some sort of simplification such as the two-stream approximation.
In searching this article also appeared, https://pdfs.semanticscholar.org/6252/b880e44b1860298521e2956c2e33eb8e3434.pdf which seems relevant.

Update.
OP is asking about source functions in scattering atmospheres.
Take a look at Chandrasekhar’s book, Radiative Transfer. Here is a screenshot which discusses the form of the source function for scattering atmospheres.

A: Derived from the Radiation transfer eqn. one has the so-called time-stationary Photon Diffusion eqn. (which is the screened Poisson eqn.)
$$D \Delta I(\vec{x})- \mu_a I(\vec{x}) + S(\vec{x})=0,$$
where $S$ is the light source, $D = \frac{1}{3(\mu_a + \mu_s)}$ is the diffusion constant and $\mu_s$ the reduced scattering coefficient. This is the usual Diffusion eqn. added by an absorption term.
1.) For an infinite planar source (xy-plane) this reduces to 
$$D I''(z) = \mu_a I(z), z>0$$
with solution $I(x,y,z)= I(z) = \frac{I_0}{2\sqrt{D \mu_a}} \exp\left(-\sqrt{\frac{\mu_a}{D}} z\right)$, which coincides with Beer-Lambert law for $\mu_s=0$ with $\sqrt{3}\mu_a$ as effective attenuation.
2.) For a line source (z-axis) this reduces to (using cylinder coordinates and symmetry)
$$D \left[I''(r) + \frac{1}{r} I'(r)\right] = \mu_a I(r), r>0$$
which is the Bessel differential eqn. with solution $I(x,y,z) = \frac{I_0}{2\pi D} K_0\left(\sqrt{\frac{\mu_a}{D}}\|(x,y)\|\right)$, where $K_0$ is the modified Bessel function of the 2nd kind of order zero. This gives the following intensity profile.

The singularity in $(x,y) = o$ makes sense. Infinite intensity in a point does not mean an 'infinite amount' of light similar as for a probability distribution (like Beta distr.). The intensity of a single light ray (without any effects) is given by a Dirac delta. 
Overall, including the attenuation of the initial ray itself gives the final result
$$I(x,y,z) = \frac{3\sqrt{3} I_0 (\mu_a + \mu_s)^{3/2}}{4\pi \sqrt{\mu_a}}\exp\left(-\sqrt{3\mu_a(\mu_a + \mu_s)} z \right) K_0\left(\sqrt{3\mu_a(\mu_a + \mu_s)}\|(x,y)\|\right).$$
Intensity in $xz$-plane.

The reduced scattering coefficient is given by $\mu_s=\mu(1-g)$, where
$$g = \frac{1}{4\pi}\int_0^{2\pi} \int_0^{\pi} |\cos\theta |(1+\cos^2 \theta)\sin \theta  d \theta d \phi = \frac{3}{4} $$
is the Rayleigh scattering anisotropy.
