## Abstract

A Monte Carlo algorithm has been developed to investigate the effects of multiple scattering on the volume scattering function measured by the LISST-VSF instrument. The developed algorithm is compared to experimental results obtained from bench-top measurements using ${508}\,\textrm{nm}$ spherical polystyrene beads and Arizona test dust as scattering agents. The Monte Carlo simulation predicts measured volume scattering functions at all concentrations. We demonstrate that multiple scattered light can be a major contributor to the detected signal, resulting in errors in the measured volume scattering function and its derived inherent optical properties. We find a relative error of 10% in the scattering coefficient for optical depths ∼0.4, and it can reach 100% at optical depths ∼2.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Multiple scattering may be a significant source of systematic errors when measuring the inherent optical properties (IOPs) of particle-rich natural waters [1,2]. IOPs include the absorption coefficient and the volume scattering function (VSF) with its derived properties, such as the scattering coefficient $b$, and the backscattering coefficient $b_b$, all of which may vary spectrally. Turbid particle-rich waters are characteristically found in coastal and estuarine regions, which are of increased interest to the optical oceanography community due to the ecological and economic importance of these regions [3–7]. Increased knowledge of IOPs in such optically complex waters would be beneficial for remote sensing observations, environmental monitoring and underwater optical communication and visibility studies. The assumption of single scattering, where each photon is only scattered once, is an important approximation when measuring the VSF or other IOPs, as this negates the use of extensive radiative transfer calculations. However, this assumption will no longer be valid if the amount of multiple scattering is too high. In a study by Chami et al. [8], it was found through radiative transfer modeling that multiple scattering may contribute to as much as $\sim 94\%$ of the radiance reflectance when the ratio of backscattering to absorption is larger than 0.3.

The VSF $\beta (\theta )$ is defined as the radiant intensity $dI$ scattered per elemental volume $dV$ in the direction $\theta$ per unit incident irradiance $E$,

Here, no azimuthal dependency is assumed. From this expression one can derive the scattering coefficient,

The scattering phase function is given as

The VSF has been sparingly measured *in situ*, especially at large scattering angles, due to lack of available instrumentation. The LISST-VSF (Sequoia Scientific) is a recently developed commercial instrument for measuring the VSF *in situ* or in the laboratory for angles 0.1-150$^{\circ }$ at a ${515}\,\textrm{nm}$ wavelength [9]. The instrument uses a ring detector for measuring the light scattered from the incident laser beam (similar to other LISST-instruments) up to a scattering angle of 14.4$^{\circ }$, while a rotating eyeball detector is used for angles 15-150$^{\circ }$. While the geometries of the two detectors are different, both make the assumption that all detected light has been scattered solely from the beam. Light lost due to absorption or secondary scattering along the path of the scattered light is corrected for (see Section 2.2.2), but no additional light is assumed to be scattered into the detectors. It has been shown that applying an absolute calibration of the eyeball detector, using a method first presented by Hu et al. [10], yields a distinct discontinuity between the ring detector VSF and eyeball detector VSF in turbid waters [1]. In addition, the scattering values seem to be unreasonably elevated compared to parallel LISST-200X measurements, for which the sample chamber has a much shorter optical path length of ${2.5}\,\textrm{cm}$, compared to the ${15}\,\textrm{cm}$ in the LISST-VSF. Microscopic polystyrene beads of precise and accurate size and concentration enable direct comparison of Mie theory and LISST-VSF measurements. Thus, polystyrene beads has been used to calibrate the LISST-VSF instrument [10–13]. In earlier laboratory measurements with 190 and ${508}\,\textrm{nm}$ diameter polystyrene beads, we have seen similar VSF discontinuities at larger bead concentrations, as well as a non-linear relationship between scattering and particle concentration for 25 $\mu$m beads [1]. The discrepancies seen in both field and laboratory measurements motivate a deeper inquiry into the multiple scattering effects on the LISST-VSF measurements.

The single scattering approximation greatly simplifies the computation of IOPs from *in situ* measurements. The optical depth $\tau = cL$ has previously been used to state whether the single scattering assumption is a good approximation or not. Here, $c$ is the attenuation coefficient, and $L$ is the distance travelled in the medium. In the case of the LISST-VSF, this distance is $L={0.15}\textrm{m}$, which is the characteristic path length of the instrument. The closely related scaled optical depth $\tau ^* = cL(1-g) = \tau (1-g)$, is also used. Here, $g$ is known as the asymmetry factor, calculated as the mean cosine of the scattering angle of the phase function [14]. In a paper by Koestner et al. [12], the general condition for single scattering is first given as $\tau << 1$. They also state that $\tau ^* << 1$ can be used to determine the regime where the single scattering approximation is valid. However, these two conditions are quite different, as the asymmetry factor g is typically between 0.89 and 0.95 for natural waters [14]. In a study by Hu et al. [10], $c = 5$ m$^{-1}$ is stated as the upper limit for single scattering, which translates to $\tau =0.75$ for the LISST-VSF instrument. By contrast, van de Hulst [15] state a significantly stricter general condition for single-scattering: For an optical depth $\tau < 0.1$ (corresponding approximately to $\tau ^* < 0.01$), single scattering can be assumed. For $0.1 < \tau < 0.3$, double-scattering corrections may be necessary, and for $\tau > 0.3$ multiple scattering must be taken into account. In a study by Agrawal and Mengüç [16], unnormalized VSF measurements were made using a nephelometer (with a fixed path length) over a wide range of mono- and polydispersed polystyrene bead concentrations. Here, the single-scattering condition was shown to agree well with the measurements up to an optical depth of 0.1, and by including an analytical double-scattering term, there was good agreement between theoretical results and measurements up to $\tau = 1.0$. Similar measurement series with increasing bead concentrations were performed by Chae and Lee [17], where a non-linear trend could be seen when scattering measurements were compared with optical thickness for $\tau >0.1$, although the authors concluded that the relationship was reasonably linear.

Monte Carlo (MC) simulations offer an excellent tool for investigating optical instrumentation due to the straight-forward implementation of geometry in the simulations. Analytical solutions to light transport can be difficult, if not impossible, to find without the use of approximations. MC simulations, are on the other hand numerical solutions and can usually be implemented without approximations. MC simulations have previously been used to investigate a variety of optical instruments. McKee et al. [18] applied a MC simulation in order to investigate the scattering collection performance of the AC-9 (WET Labs, now Sea-Bird Scientific) reflecting tube absorption meter. Based on the MC simulation, they developed an iterative scattering error correction procedure, which was later improved upon [19]. In Doxaran et al. [20] and Vadakke-Chanat et al. [21], the uncertainty and effects of absorption was analyzed and quantified for back scattered light in the ECO-BB9 (WET Labs, now Sea-Bird Scientific) sensor. Other examples can be found in [22–24]. In this paper, we present a MC simulation to investigate the effects of multiple scattering in LISST-VSF measurements. The aim of the work is to explain the observed features in the measured VSF, and to quantify errors in the phase function and scattering coefficient originating from multiple scattering.

## 2. Methods

#### 2.1 LISST-VSF measurements

Two different scattering agents are used in this study: Polystyrene beads and Arizona test dust. The polystyrene beads solution is monodisperse with a particle size distribution centered at ${508}\textrm{nm}$ and with a full-width half-maximum of ${16}\,\textrm{nm}$. The Arizona test dust solution was made from dry powder with original particle sizes ranging from 0-${50}\,\mathrm{\mu}\textrm{m}$. The mixed solution was left to stabilize overnight so that larger particles would settle, leaving only the smaller particles suspended in the solution. Consequently, the particle size distribution for the Arizona test dust measurements is not known but is expected to be dominated by the smaller particles.

All VSF measurements were performed in the laboratory using the LISST-VSF instrument in its benchtop mode. Each measurement series commenced with ${1630}\,\textrm{mL}$ of milli-Q water being poured into the instrument sample chamber. After waiting one hour for bubbles to dissipate, a blank measurement was made. Subsequently, the scattering agent was added from a master solution to the sample chamber in predetermined amounts using pipettes for precise concentrations. The milli-Q water used in the blank was also used in the following measurements, yielding an effective subtraction of pure water scattering and other possible optical losses. A series of measurements were made with subsequent additions of scattering agent solution, going from low to very high concentrations. The concentrations used in this study are given in Table 1 (Section 3.1). To minimize uncertainties connected to the particle concentrations, the VSF measurements were scaled by the ratio of theoretical (Mie theory) and measured (LISST-VSF) attenuation,

#### 2.2 Monte Carlo simulation

Here, we present a Monte Carlo simulation based on the random sampling of variables. The variables are expressed as probability distributions, so that the value of a certain variable can be sampled by generating a random number $\xi \in [\begin{matrix}{0},&{1}\end{matrix}]$. The geometry of the simulation is presented in Fig. 1. The simulation boundary is defined by a cylinder with a length stretching from $z={0}\,\textrm{cm}$ to $z={15}\,\textrm{cm}$, and a radius of ${2}\,\textrm{cm}$. The light source is positioned at the center of the top end of the cylinder ($z={15}\,\textrm{cm}$) and the photons are emitted along the z-axis towards the bottom end of the cylinder. The eyeball detector is positioned at the cylinder wall at $z={5}\,\textrm{cm}$, so that the distance from the light source is ${10}\,\textrm{cm}$ along the z-axis. The ring detectors are positioned at $z={0}\,\textrm{cm}$ and covers the entire bottom area of the cylinder. The cylinder wall act as a perfect absorber, so that all photons crossing the cylinder boundaries are immediately eliminated. It is desirable to minimize the computation time of the simulation without compromising the precision. Thus, the cylinder radius of ${2}\,\textrm{cm}$ is chosen, as it is considered statistically highly unlikely for a photon that has crossed this cylinder wall boundary to be scattered back into the sample volume and into a detector.

The individual photon’s trajectory through the sample volume is calculated following the steps outlined in the flow chart presented in Fig. 2. The distance between photon-particle interaction will here be referred to as step size, where the step size distribution can be derived from the attenuation of a beam,

Here, $I_0$ is the initial intensity of the beam, $c$ is the attenuation coefficient, $l$ is the distance over which the beam is attenuated, and $I$ is the intensity after traveling the distance $l$. A random step size can then be sampled from the distribution where the sampled step size $s$ can be calculated as

The new position of the photon is calculated based on the traveling direction and the sampled step size, and the scattering angle distribution is obtained from the scattering probability density function (PDF) of the sample in question. The scattering PDF is given as

and integration yields where $p(\theta )$ is the scattering phase function. The scattering angle $\theta$ is randomly sampled by generating a new $\xi \in [\begin{matrix}{0},&{1}\end{matrix}]$. The sampled scattering angle $\theta _{sc}$ is then found from the relationA lookup table containing the values of $\theta _{sc}$ for different $\xi$ is made prior to the simulation. The azimuthal angle $\phi$ is sampled from a uniform distribution from $0$ to $2\pi$. When performing a simulation for the ${508}\,\textrm{nm}$ spherical polystyrene beads, a theoretical phase function computed from Mie theory is used as input. For Arizona test dust, the input phase function is obtained from a LISST-VSF measurement at low concentration, so that one can assume that the single scattering approximation is valid.

In order to account for absorption by the sample volume, we give the photons an initial weight $w=1$, and at each scattering event, the weight is reduced according to

where $w_{i+1}$ is the weight after scattering, $w_{i}$ is the weight before scattering, $c$ is the attenuation coefficient, and $b$ is the scattering coefficient. This avoids having to initiate a new photon every time one is absorbed. As the LISST-VSF measurement is calibrated so that the absorption and scattering from the water itself is zero, there is no need to implement the attenuation from water in the simulation.### 2.2.1 Simulation of detectors

### 2.2.1.1 Ring detectors

In order to detect the scattered intensity in the range $0.09-{14.4}^{\circ}$, the LISST-VSF uses a set of ring detectors, covering different parts of the range. A lens is positioned above the ring detectors, focusing the light at the different rings depending on the angle of incidence. This detection scheme is simulated by simply defining a flat circle surface, corresponding to the lens, and registering the angle of incidence onto this surface, as illustrated in Fig. 1(a) . The registered photons are then binned according to their angle of incidence, matching the bin sizes of the LISST-VSF.

### 2.2.1.2 Eyeball detector

The eyeball detection scheme is illustrated in Fig. 1(b). The light passes through an opening in the eyeball, leading to a set of reflecting mirrors. The mirrors guide the light to a spatial filter consisting of two lenses and a pinhole, resulting in an acceptance angle of $2\theta _a={0.9}^{\circ}$. The total power reaching the detector is thus dependent on the opening of the eyeball and the acceptance angle. While the LISST-VSF eyeball rotates in order to measure the angle dependency of the scattered light, the simulation can detect scattered light from all angles $\theta$ simultaneously. This is achieved by defining a sphere positioned at the center of the eyeball, where passing through this sphere is the equivalent of passing through the eyeball opening and being within the acceptance angle of the spatial filter. The principle behind this design is illustrated in Fig. 1(c), and while this is a major simplification of the eyeball detection system, the results presented in this paper demonstrate that it works very well. A sphere radius of ${3}\,\textrm{mm}$ was found to give a good fit to the experimental results. The simulation also allows us to extend this eyeball design around the cylinder, creating a ring torus (see Fig. 1(a)). This is the equivalent of having eyeball detectors covering the entire circumference of the cylinder, as opposed to the LISST-VSF which only has one eyeball detector at a single location. Thus, the simulation can detect scattered photons for all azimuth angles.

### 2.2.2 Processing simulated data

To obtain a VSF from the MC simulation, the raw data, i.e. number of photons detected at each angle, must be processed. The VSF is calculated based on the method presented in [25]. The number of scattered photons detected by the ring detectors can be calculated for each bin (ring) as

where $e^{-cl}$ is the attenuation factor and $l$ is the total distance traveled from the laser injection point to the detector, assuming single scattering. The factor $x$ is the length of the laser beam contributing to the signal (see Fig. 3), $N_0$ is the total number of simulated photons, equivalent to the laser power, $\beta _i(\theta )$ is the VSF for bin $i$, and $\theta _{i,h}$ and $\theta _{i,l}$ is the high and low limit of the bin, respectively. The symbol $\phi$ does here refer to the fraction of a full circle covered by the ring detectors. For practical reasons, the ring detectors in the LISST-VSF only cover $1/6$ of the full circle, meaning that $\phi =1/6$. In the simulation however, we can use the full circle so that $\phi =1$. Rearranging Eq. (11), we get the following expression for the VSFHere, the factor $e^{cl}$ accounts for photons lost along the path to the detector. For the ring detector, $l$ does not vary much with detection angle, and is set to $l=L={0.15}\,\textrm{m}$. An identical correction for attenuation is done for the LISST-VSF measurements. For angles smaller than $\arctan (r/L)$, where $r$ is the radius of the lens focusing light at the ring detectors, the entire beam contributes to the signal, so that $x=L={0.15}\,\textrm{m}$. For angles larger than this, only a fraction of the beam contributes as the photons scattered from the top end of the beam does not hit the lens (see Fig. 3(a)). In this case, $x$ is calculated as $x=r/\tan {\theta _i}$, such that Eq. (12) becomes

The same approach is used to calculate the VSF from the signal obtained by the eyeball detector. In this case, only a small fraction of the beam contributes to the signal, where $x$ is calculated as $x=d/\sin {\theta _i}$ (see Fig. 3(b)). Here, $d={6}\,\textrm{mm}$ is the diameter of the eyeball detection sphere defined in Section 2.2.1.2. The VSF for the eyeball detector can then be calculated for bin $i$ as,

Here, the path length $l$ is dependent on the angle of detection, and can be calculated as

where $l_{\textrm {eye}}={0.10}\,\textrm{m}$ is the distance from the top end of the cylinder to the eyeball detector along the z-axis, and $R$ is the distance from the laser beam to the eyeball detector measured at a right angle. Again, we can set $\phi =1$ due to the full circle coverage of the eye detector.In order to ensure a similar signal for all samples, the simulation is set to run until $10^7$ photons are detected by the eye detector, i.e. $\sum N_{i, eye}=10^7$. The same limit is set for the ring detector, at which the detector stops counting detected photons. The final signal in the ring detector is then calculated as $N_i=N \cdot N_{i,\textrm {stop}}/N_{\textrm {stop}}$, where $N_{i}$ is the adjusted number of photons in bin $i$, $N$ is the total number of photons simulated, $N_{i,\textrm {stop}}$ is the number of photons in bin $i$ at which the limit of $10^7$ detected photons is reached, and $N_{\textrm {stop}}$ is the total number of simulated photons when the number of detected photons reaches $10^7$.

## 3. Results and discussion

#### 3.1 Monte Carlo simulations compared with experimental results

The results obtained using the polystyrene beads as the scattering agent are presented in Fig. 4. The simulation input parameters $b$ and $c$ are given in Table 1, along with the corresponding experimental particle concentrations $C$ and optical depths $\tau$. Measurements by the LISST-VSF are compared with single-scattering theoretical values computed from Mie theory, and simulation results from the Monte Carlo algorithm. At low concentrations, the theoretical, measured, and simulated VSF are in overall good agreement. For increasing particle concentrations, the measured VSF increasingly deviates from the theoretical VSF, both in magnitude and shape. For the highest bead concentration, the amount of measured scattering at large angles is more than an order of magnitude larger than predicted from single scattering, and the oscillating features of the ${508}\,\textrm{nm}$ bead phase function have vanished. Another prominent feature at higher concentrations is the discontinuity at ${15}^{\circ}$. This angle is the boundary between where the VSF is measured by the ring detector and where the eye detector is used.

By contrast, the simulation predicts a VSF with closer resemblance to the LISST-VSF measurements, even at the highest scattering coefficient $b\approx 33$ m$^{-1}$. Some minor deviations can be seen. At small scattering angles ($<{7}^{\circ}$) and low concentrations, the LISST-VSF measured enhanced forward scattering relative to theory, whereas the simulated and theoretical curves are almost indistinguishable. The reason is not entirely clear, but is likely attributed to bead flocs, i.e., beads clumped together and effectively acting as larger particles [1,26]. Increased forward scattering for larger particles is consistent with Mie theory. The bead flocs may then deflocculate (break apart) when mixing with water during the experiments, resulting in the spike gradually vanishing. The concentrated polystyrene bead solution was also shaken before each addition of beads to the LISST-VSF sample chamber, so later additions of beads were more thoroughly shaken. Thus, the later additions of beads may have contained less flocculated beads, adding to the reduction of the forward spike. There are other possible explanations for the observed spike in forward scattered intensity. Local heating of the sample volume in the vicinity of the laser beam might cause a temperature gradient to arise, which in turn results in density gradients and deviations in the refractive index [27]. Such deviations can cause scattering of light, and thereby an increase in the measured forward scattering. However, if this were the case, one would also expect to measure increased attenuation, which is not observed in our experiments. Thus, this explanation is considered unlikely. The discontinuity at $15^{\circ }$ is also apparent in the simulated VSF, but appears larger in the experimental data for sample Nos. 6-9. After inspection of the experimental raw data, we find that this is due to saturation of the eye detector sensor at the smallest angles ($15-{25}^{\circ}$). At the highest concentration, the saturation is quenched by increased attenuation of the light (which is compensated for in the data processing), hence the improved agreement between the experiment and simulation for sample No. 10. From 140° and onward, the simulated backscattering increases relative to the experimental backscattering, with increasing particle concentration. This is most likely due to the simulated eyeball detection being a simplified version of the instrument eyeball detection. Consequently, the two different eyeball detection schemes might probe the sample volume slightly different, resulting in deviations in the VSF when the turbidity becomes very high.

The Arizona test dust results are presented in Fig. 5. The simulation input parameters $b$ and $c$ are given in Table 1, with the corresponding experimental relative particle concentrations (specific particle concentration is not known) and optical depths $\tau$. Here, the theoretical VSF is acquired from a low concentration sample and scaled according to the relative concentrations $C_r$ presented in Table 1. The results follow the same trend as for the polystyrene beads, showing clear deviations between the predicted values for single scattering and the measurements at high concentrations, and close agreement between simulated and experimental results. A clear difference in the two phase functions, are the oscillatory features seen in the ${508}\,\textrm{nm}$ beads phase function (see Fig. 4), as opposed to Arizona test dust phase function, where such features are not present. The oscillatory features appear due to constructive and destructive interference of the electromagnetic field scattered by individual particles [12,28]. The angular positions of constructive and destructive interference are dependent on particle size and refractive index, such that the monodisperse ${508}\,\textrm{nm}$ beads samples will display oscillations in the phase function, while the polydisperse Arizona test dust samples will not.

Results seen in Figs. 4 and 5 validate the developed Monte Carlo simulation of the LISST-VSF, and offer robust explanations of instrument artifacts seen in both laboratory and field measurements. Multiple scattering causes both erroneously large detected scattering and an altered phase function. The discontinuity at ${15}^{\circ}$ arises due to the sudden decrease in path length going from ring detectors to eyeball detector (see Fig. 1). This decrease in path length has a double effect on the measured VSF. Firstly, longer path lengths result in more multiple scattered photons. Secondly, the attenuation factor $e^{cl}$ is larger for the ring detectors (Eq. (12)) than for the eye detector (Eq. (14)) in the forward direction. Thus, the increase in number of detected photons $N_i$ due to multiple scattering contributes more to the absolute error in the VSF measured by the ring detector compared to the one measured by the eyeball detector. However, this does not impact the relative error. Saturation of the eyeball detector may also be a significant factor in waters dominated by scattering processes.

#### 3.2 Error analysis

### 3.2.1 Convergence test

A convergence test has been performed for the ${508}\,\textrm{nm}$ beads sample No. 1 (see Table 1 for details about the sample). Convergence is tested for the scattering coefficient $b$, calculated from Eq. (2), and the results are presented in Fig. 6. For the test, the simulation is set to run until a specific number of photons $N_{\textrm{eye}}$ are detected by the eye detector, starting at 1221 photons and doubling for each simulation up til $10^7$ (same as for the results presented in Figs. 4 and 5). The simulation is repeated 10 times for each $N_{\textrm{eye}}$, providing enough data to calculate a coefficient of variance (CV), also known as relative standard deviation. The CV is calculated as $c_v=\frac {\sigma }{\mu }\times 100\%$, where $\sigma$ is the standard deviation of the mean $\mu$. The calculated scattering coefficients $b$ from the simulations are presented in Fig. 6(a), and the CV of $b$ is plotted in Fig. 6(b), in addition to a theoretical fit calculated as $c_{v}=C/\sqrt {N_{\textrm{eye}}}$, where $C$ is a proportionality constant.

From Fig. 6, one can see that $b$ converges relatively fast. In fact, the CV is below $1\%$ when $N_{\textrm{eye}}$ is larger than $10^4$. While the simulated VSF obtained from a signal of $10^4$ detected photons would be very noisy, the scattering coefficient can be accurately predicted as it is obtained from integrating the VSF over the angels ${0.09}^{\circ}-{150}^{\circ}$ (see Eq. (2)). At $10^7$ detected photons, the CV is approximately 0.04%, giving a very high precision. A similar convergence test was performed for sample number 5, and no significant difference in CV was found. Thus, it was concluded that the optical depth has little to no influence on the CV for the scattering coefficient.

### 3.2.2 Multiple scattering error

The validity of the single-scattering approximation may readily be assessed by comparing the Monte Carlo simulation results with the predicted results from assuming single-scattering. For the error analysis, the theoretical scattering coefficients and phase functions are calculated from the VSFs labeled "Theory" in Fig. 4 and 5. The theoretical, experimental, and simulated scattering coefficients, denoted $b_t$, $b_e$, and $b_s$, respectively, are given in Table 2. The scattering coefficients are calculated by integrating from ${0.09}^{\circ}$ to ${150}^{\circ}$, so that theory, experiment and simulation are considered over the same angular spectrum. Thus, the theoretical scattering coefficient is not the same as the input scattering coefficient given in Table 1.

From Table 2, one can see that a difference in the theoretical single-scattering scattering coefficient $b_t$ and the experimental scattering coefficient $b_e$ can be quite dramatic for high particle concentrations. As the LISST-VSF attenuation measurements have been shown to have a high accuracy and are virtually unaffected by multiple scattering [1,23], one would see a similar dramatic difference in the measured absorption $a$, calculated as $a=c-b$. In fact, in some cases, when $b>c$, one would get a negative absorption coefficient. The relative errors due to multiple scattering for the scattering coefficient and phase function are plotted in Fig. 7. Three different relative errors are presented. The desired output from a LISST-VSF measurement is the single-scattering phase function and scattering coefficient, and are thus taken as true values when computing the error in the LISST-VSF measurements. This is also the case when comparing simulated results to single scattering. When comparing simulation and experiment, the experimental results are taken as true values, as the goal of the simulation is to reproduce the LISST-VSF measurements. The scattering coefficient relative error $e_b$ is computed from

For the phase function, where the error may vary significantly with angle, the mean relative error $e_{p}$ is calculated, using

The measured phase function is measured at log-spaced angles at $<{15}^{\circ}$. To avoid bias towards forward scattering effects, the data was interpolated to evenly spaced angles.

The scattering coefficient and phase function errors for the ${508}\,\textrm{nm}$ beads samples are plotted in Figs. 7(a) and (b), respectively, and the errors for the Arizona test dust samples are plotted in Figs. 7(c) and (d). The error is plotted against the dimensionless parameter $\tau$. As seen in Fig. 7, the simulation predicts a scattering coefficient and phase function close to those derived from the LISST-VSF measurements. For the scattering coefficients, the relative error is always $<13\%$, and $<10\%$ for the phase function. As discussed in section 3.1, a large portion of the error in the simulation relative to the LISST-VSF measurements can be attributed to measurement errors in the experimental results. Comparing simulation and experiment to the single-scattering theoretical values, the errors in the scattering coefficients arise solely due to the added intensity from multiple scattered photons. The phase function, however, is normalized so that the addition of multiple scattered photons does not alone cause the error. Rather, it is the angular distribution of multiple scattered photons that causes the error. Due to the phase functions being highly forward peaked, the photons scattered more than once tend to be distributed to angles close to the original scattering angle. As a consequence, phase functions that varies rapidly with scattering angle is more subjected to multiple scattering errors in the phase function, as the relative increase in signal becomes large when a low signal part of the spectrum neighbors a high signal part. This is especially evident looking at the relative error in the phase functions (Figs. 7(b) and (d)), where the error is significantly larger for the ${508}\,\textrm{nm}$ beads than for the Arizona test dust. This is due to the low signal dips in the ${508}\,\textrm{nm}$ beads phase function, approximately at ${60}^{\circ}$ and ${110}^{\circ}$, gaining a relatively large quantity of multiple scattered photons from the neighboring higher signal parts. Ultimately, as the optical depth approaches infinity, the angular distribution of photons becomes completely random, resulting in a uniform VSF.

Comparing simulation to theory, the error in the scattering coefficient is calculated to be $e_b=0.4\%$ for the ${508}\,\textrm{nm}$ at an optical depth of $\tau =0.027$. The error increases with optical depth reaching $e_b=29\%$ for $\tau =1.1$, and finally $e_b=278\%$ for $\tau =4.9$. As the CV was calculated to be $c_v=0.04\%$, the numerical uncertainty in the simulated results can be considered negligible, except for at the smallest optical depths. The error in the phase function is, for the most part, significantly smaller than for the scattering coefficient. The exceptions are for the three smallest optical depths. This is probably due noise being the dominating source of error when the difference between simulated and theoretical VSF becomes very small. The scattering coefficient is not subjected to noise in the same way, as the scattering coefficient is calculated by integrating the VSF over the entire angular range (see Eq. (2)). At the smallest optical depth $\tau =0.027$, the error is calculated to be $e_p=2\%$, increasing to $e_p=123\%$ for the largest optical depth of $\tau =4.9$. For the Arizona test dust, the error in the scattering coefficient goes from $e_b=2\%$ to $e_b=250\%$ when the optical depth goes from $\tau =0.013$ to $\tau =3.9$. The error in the phase function goes from $e_p=1\%$ to $e_p=40\%$ for the same values of $\tau$. An overview of the optical depths at which the scattering coefficient and phase function reach specific relative errors is presented in Table 3. The values presented are calculated by linear interpolation of the results labeled as "Simulation vs. theory" in Fig. 7.

As seen in the results presented here, a ${15}\,\textrm{cm}$ path length is sufficient to cause considerable errors in the measured VSF and its derived inherent properties, here represented by the scattering coefficient and phase function. Thus, it becomes interesting to compare the LISST-VSF with other optical equipment and measurements. For instance, the historical Petzold measurements are often used as a reference to *in situ* data [29]. For these measurements, two different instruments were used, depending on the scattering angle. The instrument used to measure the low angle scattering (approximately ${0.06}^{\circ}-{0.4}^{\circ}$) had a path length of ${50}\,\textrm{cm}$, while the range ${10}^{\circ}-{170}^{\circ}$ was measured at a path length of ${18.8}\,\textrm{cm}$. The VSF was measured for a range of different ocean waters, with an attenuation coefficient peaking at 2.190, which translates to an optical depth of 1.1, for the longest path length of ${50}\,\textrm{cm}$. From Table 3, one can see that this results in a relative error of approximately 30% in the measured scattering coefficient, assuming the equipment is subjected to multiple scattering similarly to the LISST-VSF.

## 4. Summary and conclusion

We have demonstrated that the errors originating from multiple scattering in the LISST-VSF measurements can be significant for large optical depths, which largely explains *in situ* LISST-VSF measurement artifacts seen in turbid waters measurements. The developed Monte Carlo algorithm has been validated by employing ${508}\,\textrm{nm}$ polystyrene beads and Arizona test dust as scattering agents, and has proven to be an accurate tool for analyzing errors in both the measured phase function and scattering coefficient. A convergence test was performed for the MC algorithm, from which we found a coefficient of variance of $c_v\leq 0.04\%$ for the scattering coefficient. Comparing simulation to experiment, the relative error is <13% for the scattering coefficient and <10% for the phase function, for the concentrations and scattering agents tested in this study. A large fraction of these errors can be attributed to errors in the LISST-VSF measurements, caused by particle flocculation, or saturation of the eyeball detector.

Comparing the simulated VSF to the single-scattering VSF, the relative error in the scattering coefficient reaches $10\%$ for an optical depth of $\tau \approx 0.4$ for both scattering agents, while an error of $100\%$ is reached for $\tau \approx 2.8$ and $\tau \approx 2.0$ for the ${508}\,\textrm{nm}$ beads and Arizona test dust, respectively. The error in the phase function was found to be significantly larger for the ${508}\,\textrm{nm}$ beads than for the Arizona test dust. A $10\%$ error is found for $\tau \approx 0.48$ for the ${508}\,\textrm{nm}$ beads, reaching $100\%$ for $\tau \approx 4.1$. The Arizona test dust, however, does not reach a $10\%$ error before $\tau \approx 1.2$ and does not go above $50\%$ for the concentrations investigated. The large difference in error between the two phase functions is attributed to the finer details in the ${508}\,\textrm{nm}$ beads phase function. Thus, we have shown that errors in both the scattering coefficient and phase function are dependent on the phase function, but to what extent is still not fully explored. Furthermore, both sample sets investigated in this paper are heavily dominated by scattering, increasing the ratio of absorption to scattering is expected to influence the error when plotted against optical depth.

## Acknowledgments

We would like to thank one of our anonymous reviewers who brought to our attention the possibility of multiple scattering errors in the Petzold measurements.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **H. Sandven, A. S. Kristoffersen, Y.-C. Chen, and B. Hamre, “In situ measurements of the volume scattering function with LISST-VSF and LISST-200X in extreme environments: evaluation of instrument calibration and validity,” Opt. Express **28**(25), 37373–37396 (2020). [CrossRef]

**2. **Y. C. Agrawal and H. C. Pottsmith, “Instruments for particle size and settling velocity observations in sediment transport,” Mar. Geol. **168**(1-4), 89–114 (2000). [CrossRef]

**3. **C. B. Mouw, S. Greb, D. Aurin, P. M. DiGiacomo, Z. Lee, M. Twardowski, C. Binding, C. Hu, R. Ma, T. Moore, W. Moses, and S. E. Craig, “Aquatic color radiometry remote sensing of coastal and inland waters: Challenges and recommendations for future satellite missions,” Remote. Sens. Environ. **160**, 15–30 (2015). [CrossRef]

**4. **M. Chami, E. B. Shybanov, G. A. Khomenko, M. E.-G. Lee, O. V. Martynov, and G. K. Korotaev, “Spectral variation of the volume scattering function measured over the full range of scattering angles in a coastal environment,” Appl. Opt. **45**(15), 3605–3619 (2006). [CrossRef]

**5. **M. Soja-Woźniak, M. Baird, T. Schroeder, Y. Qin, L. Clementson, B. Baker, D. Boadle, V. Brando, and A. D. L. Steven, “Particulate backscattering ratio as an indicator of changing particle composition in coastal waters: Observations from Great Barrier Reef waters,” J. Geophys. Res.: Oceans **124**(8), 5485–5502 (2019). [CrossRef]

**6. **V. G. Harharasudhan and P. Shanmugam, “Modelling the particulate backscattering coefficients of turbid and productive coastal waters,” Ocean Sci. J. **54**(2), 147–164 (2019). [CrossRef]

**7. **C. Nima, Ø. Frette, B. Hamre, S. R. Erga, Y.-C. Chen, L. Zhao, K. Sørensen, M. Norli, K. Stamnes, and J. J. Stamnes, “Absorption properties of high-latitude Norwegian coastal water: The impact of CDOM and particulate matter,” Estuarine, Coastal Shelf Sci. **178**, 158–167 (2016). [CrossRef]

**8. **M. Chami, D. McKee, E. Leymarie, and G. Khomenko, “Influence of the angular shape of the volume-scattering function and multiple scattering on remote sensing reflectance,” Appl. Opt. **45**(36), 9210–9220 (2006). [CrossRef]

**9. **W. H. Slade, Y. C. Agrawal, and O. A. Mikkelsen, “Comparison of measured and theoretical scattering and polarization properties of narrow size range irregular sediment particles,” in 2013 OCEANS-San Diego, (IEEE, 2013), pp. 1–6.

**10. **L. Hu, X. Zhang, Y. Xiong, and M.-X. He, “Calibration of the LISST-VSF to derive the volume scattering functions in clear waters,” Opt. Express **27**(16), A1188–A1206 (2019). [CrossRef]

**11. **T. Harmel, M. Hieronymi, W. Slade, R. Röttgers, F. Roullier, and M. Chami, “Laboratory experiments for inter-comparison of three volume scattering meters to measure angular scattering properties of hydrosols,” Opt. Express **24**(2), A234–A256 (2016). [CrossRef]

**12. **D. Koestner, D. Stramski, and R. A. Reynolds, “Measurements of the volume scattering function and the degree of linear polarization of light scattered by contrasting natural assemblages of marine particles,” Appl. Sci. **8**(12), 2690 (2018). [CrossRef]

**13. **D. Koestner, D. Stramski, and R. A. Reynolds, “Polarized light scattering measurements as a means to characterize particle size and composition of natural assemblages of marine particles,” Appl. Opt. **59**(27), 8314–8334 (2020). [CrossRef]

**14. **M. Jonasz and G. Fournier, * Light Scattering by Particles in Water: Theoretical and Experimental Foundations* (Elsevier, 2007).

**15. **H. C. van de Hulst, * Light Scattering by Small Particles* (Courier Corporation, 1981).

**16. **B. M. Agrawal and M. P. Mengüç, “Forward and inverse analysis of single and multiple scattering of collimated radiation in an axisymmetric system,” Int. J. Heat Mass Transfer **34**(3), 633–647 (1991). [CrossRef]

**17. **S.-K. Chae and H. S. Lee, “Determination of radiative transport properties of particle suspensions by a single-scattering experiment,” Aerosol Sci. Technol. **18**(4), 389–402 (1993). [CrossRef]

**18. **D. McKee, J. Piskozub, and I. Brown, “Scattering error corrections for in situ absorption and attenuation measurements,” Opt. Express **16**(24), 19480–19492 (2008). [CrossRef]

**19. **D. McKee, J. Piskozub, R. Röttgers, and R. A. Reynolds, “Evaluation and improvement of an iterative scattering correction scheme for in situ absorption and attenuation measurements,” J. Atmospheric Ocean. Technol. **30**(7), 1527–1541 (2013). [CrossRef]

**20. **D. Doxaran, E. Leymarie, B. Nechad, A. Dogliotti, K. Ruddick, P. Gernez, and E. Knaeps, “Improved correction methods for field measurements of particulate light backscattering in turbid waters,” Opt. Express **24**(4), 3615–3637 (2016). [CrossRef]

**21. **S. Vadakke-Chanat, P. Shanmugam, and B. Sundarabalan, “Monte Carlo simulations of the backscattering measurements for associated uncertainty,” Opt. Express **26**(16), 21258–21270 (2018). [CrossRef]

**22. **D. Stramski and J. Piskozub, “Estimation of scattering error in spectrophotometric measurements of light absorption by aquatic particles from three-dimensional radiative transfer simulations,” Appl. Opt. **42**(18), 3634–3646 (2003). [CrossRef]

**23. **J. Piskozub, D. Stramski, E. Terrill, and W. K. Melville, “Influence of forward and multiple light scatter on the measurement of beam attenuation in highly scattering marine environments,” Appl. Opt. **43**(24), 4723–4731 (2004). [CrossRef]

**24. **J. T. O. Kirk, “Monte Carlo modeling of the performance of a reflective tube absorption meter,” Appl. Opt. **31**(30), 6463–6468 (1992). [CrossRef]

**25. **Y. C. Agrawal, “The optical volume scattering function: Temporal and vertical variability in the water column off the New Jersey coast,” Limnol. Oceanogr. **50**(6), 1787–1794 (2005). [CrossRef]

**26. **W. H. Slade, E. Boss, and C. Russo, “Effects of particle aggregation and disaggregation on their inherent optical properties,” Opt. Express **19**(9), 7945–7959 (2011). [CrossRef]

**27. **O. A. Mikkelsen, T. G. Milligan, P. S. Hill, R. J. Chant, C. F. Jago, S. E. Jones, V. Krivtsov, and G. Mitchelson-Jacob, “The influence of schlieren on in situ optical measurements used for particle characterization,” Limnol. Oceanogr.: Methods **6**(3), 133–143 (2008). [CrossRef]

**28. **M. J. Berg, C. M. Sorensen, and A. Chakrabarti, “Explanation of the patterns in Mie theory,” J. Quant. Spectrosc. Radiat. Transfer **111**(5), 782–794 (2010). [CrossRef]

**29. **T. J. Petzold, “Volume scattering functions for selected ocean waters,” Tech. rep., Scripps Institution of Oceanography La Jolla Ca Visibility Lab (1972).